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Abstract 


The  primary  objective  of  this  research  was  to  develop  and  apply  large  eddy  simu¬ 
lation  (LES)  technology  to  some  urgent  flow  and  heat  transfer  problems  in  propulsion 
systems  and  to  contribute  to  the  physical  understanding  of  such  flows.  The  work  was 
motivated  by  the  observation  that  the  design  goals  of  high  specific  power  and  thrust 
and  low  specific  fuel  consumption  have  been  reached,  in  part,  by  an  increase  in  turbine 
inlet  temperature  and  future  improvements  in  engine  efficiency  will  place  even  greater 
demands  on  blade  cooling  procedures.  However,  current  design  codes  are  somewhat 
limited  in  accuracy  due  to  uncertainty  associated  with  modeling  for  turbulent  flow.  The 
research  was  concerned  with  both  the  film  cooling  of  external  blade  surfaces  and  the 
complex  flows  in  internal  cooling  passages.  Studies  have  been  completed  showing  the 
effects  of  rotation  on  the  heat  transfer  and  flow  in  smooth  and  ribbed  channels  and  in 
a  duct  of  square  cross-section.  A  scheme  was  developed  for  including  the  effects  of 
freestream  turbulence  on  boundary  layer  development.  Preliminary  LES  results  have 
been  obtained  for  a  single  hole  film  cooling  configuration. 


1  Introduction 


This  report  summarizes  the  research  accomplishments  achieved  under  Grant  AFOSR  F49620- 
01-1-0113,  “Application  of  Large  Eddy  Simulation  to  Cooling  and  Flow  Problems  in  Aeropropul- 
sion  Systems.”  The  objective  of  the  research  was  to  expand  the  capabilities  of  LES  technology  to 
contribute  to  the  solution  of  urgent  problems  in  propulsion  systems  and  to  contribute  to  the  physical 
understanding  of  such  flows.  The  long-term  goal  of  the  thrust  is  to  advance  large  eddy  simulation 
(LES)  technology  to  the  point  where  it  is  effective  for  computing  realistic  flows  in  a  rotating,  mul¬ 
tistage  turbomachine.  First,  however,  some  building  block  simulations  need  to  be  carried  out  for  a 
single  blade  element. 

Flows  in  propulsion  systems  are  accompanied  by  effects  such  as  rotation  and  high  turbulence 
levels  that  alter  the  basic  structure  of  turbulent  flows  in  ways  that  are  not  accurately  represented 
by  current  design  codes  that  rely  heavily  on  turbulence  models.  Significant  improvements  in  ef¬ 
ficiency  can  be  achieved  if  the  uncertainty  associated  with  design  predictions  can  be  reduced.  As 
computer  technology  continues  to  advance,  it  is  anticipated  that  large  eddy  simulation  (LES)  in 
which  turbulence  modeling  plays  a  minor  part  will  play  an  ever-increasing  role  in  the  understanding 
of  such  complex  flows  and  will  enable  preliminary  design  decisions  to  be  based  on  more  accurate 
computational  models.  Turbine  cooling  issues  have  become  very  important  since  the  design  goals 
of  high  specific  power  and  thrust  and  low  specific  fuel  consumption  have  been  reached,  in  part,  by 
an  increase  in  turbine  inlet  temperature.  Future  improvements  in  engine  efficiency  will  place  even 
greater  demands  on  blade  cooling  procedures.  Accordingly,  the  the  research  has  focused  on  devel¬ 
oping  and  applying  large  eddy  simulation  (LES)  technology  to  idealized  configurations  for  both  the 
film  cooling  of  external  blade  surfaces  and  the  complex  flows  in  internal  cooling  passages. 

Large  eddy  simulation  uses  large  computer  resources  but  is  capable  of  achieving  very  realistic 
results  because  very  little  ad  hoc  modeling  is  employed  to  represent  the  effects  of  turbulence.  The 
unsteady,  three-dimensional  motion  of  the  large  eddies  is  resolved  from  first  principles  and  model¬ 
ing  is  only  used  to  account  for  the  effects  of  eddies  smaller  than  the  computational  grid  itself.  Such 
small  eddies  tend  to  be  nearly  isotropic  and  universal.  The  role  of  LES  in  dealing  with  propulsion 
technology  is  two  fold.  First,  the  technique  can  be  used  to  supplement  the  fundamental  information 
obtained  from  experiments  in  that  effects  not  easily  established  experimentally  such  as  large  tem- 
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perature  differences,  heat  flux  levels  and  rotation  effects  can  be  imposed  in  the  simulations.  Second, 
LES  has  the  potential  to  become  a  valuable  predictive  tool  for  practical  turbulent  flows  as  computer 
technology  continues  to  advance.  For  this  to  occur,  however,  considerable  work  is  required  to  find 
innovative  ways  of  simulating  flows  at  higher  Reynolds  numbers,  computationally  efficient  ways 
of  accommodating  geometries  of  increasing  complexity,  and  for  implementing  general  inflow  and 
outflow  conditions. 

The  numerical  formulation  used  in  the  research  is  based  on  a  finite  volume  discretization  of 
the  Favre-filtered  compressible  Navier-Stokes  equations.  Significant  levels  of  heat  transfer  occur  in 
the  gas  flows  of  interest  and  the  compressible  formulation  enables  the  effects  of  property  variations 
to  be  taken  into  account.  The  scheme  solves  for  the  primitive  variables  of  velocity,  pressure,  and 
temperature  in  a  fully  coupled  manner.  The  results  reported  here  employed  a  modified  version 
of  the  LU-SGS  implicit  scheme  [1,2]  that  has  been  optimized  for  parallel  implementation  using 
the  message-passing  interface  (MPI).  Low  Mach  number  preconditioning  [3]  was  used  to  enable 
solutions  to  be  obtained  efficiently  at  the  low  Mach  numbers  arising  in  many  applications  of  interest. 
A  dynamic  subgrid-scale  model  has  been  used  in  all  simulations  reported  here. 

Facets  of  the  research  have  included: 

•  Channel  flows  with  rotation  and  with  and  without  ribs 

•  Square  duct  flows  with  rotation 

•  Algorithm  developments  including  ways  to  develop  inflow  conditions  for  the  turbulent  bound¬ 
ary  layer 

•  The  effects  of  freestream  turbulence  on  turbulent  boundary  characteristics 

•  Film  cooling 

The  organization  of  this  report  is  as  follows.  First  the  conservation  equations  will  be  presented 
in  a  general  form  that  served  as  a  starting  point  for  development  of  the  various  specific  forms  needed 
to  address  the  several  applications  listed  above.  Next  the  application  areas  will  be  discussed  in  the 
order,  rotating  channel  flows,  internal  passage  flows  including  effects  of  rotation  and  centrifugal 
buoyancy,  turbulent  boundary  layer  flows  developing  under  the  influence  of  freestream  turbulence, 
and  film  cooling  flows.  The  additions  and  changes  needed  to  specialize  the  conservation  principles 
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for  each  application  will  be  explained  in  the  separate  sections  devoted  to  that  applications.  Likewise, 
relevant  numerical  details  will  given  as  the  application  areas  are  discussed. 
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2  Governing  Equations-General 


In  this  report,  flow  variables  are  nondimensionlized  as  follows. 


x*  u*  t*  p* 

— —  Ui  —  —  t  =  - - - - r  p  =  — 

Kef  Uref  (J-ref  /Uref)  P  ref 


P*  T  =  —  e=  —  p=-^~ 

PrefU^j  Tref  U}ef  Href 

k*  C  *  R*  C  * 

Kef  (U ref /Kef)  (P/ef/Kef)  P  P/ef/Kef) 


(1) 


where  the  dimensional  variables  have  been  denoted  with  a  superscript  asterisk.  The  reference  Mach 
number  is  Mref  =  Uref/ y/y R*Tref-  The  fluid  is  assumed  to  be  an  ideal  gas  and  the  non-dimensional 
equation  of  state  is  p  =  p RT.  The  non-dimensional  coefficients  of  viscosity  and  thermal  conductiv¬ 
ity  were  evaluated  as:  /j  =  k  =  Ta  where  a  is  assumed  to  be  0.71.  The  specific  heats  and  molecular 
Prandtl  number  were  treated  as  constants. 

In  compressible  flow,  fluid  properties  may  vary  temporally  and  spatially  due  to  the  simultane¬ 
ous  fluctuations  in  dilatation,  heat  transfer  and  momentum  transfer.  A  mass-weighted  averaging 
is  recommended  for  properties  carried  by  fluid  elements  which  are  related  to  mass  conservation. 
Favre-filtering  [21]  was  utilized  for  this  purpose.  Consequently,  the  Favre-filtered  compressible 
Navier-Stokes  equations  are 


dp  d(puj)  _ 

dt  dxj 

d(fiUj)  d(pujUj)  _  dp  defy  <hjj 

dt  dxj  dxi  dxj  dxj 

d{pE)  ]  d{(pE  +  p)Uj]  _d(uiaij)  dqj  dQj 
dt  dxj  dxj  dxj  dxj 

The  overline,  (•),  denotes  a  Favre-filtered  quantity.  E  =  cvt  +  \uiUj,  and  dy  = 
Moreover,  the  equation  of  state  becomes 


(2) 

(3) 


(4) 


Rer(Sij  3  Skkdij)- 


p  =  RpT 


(5) 


The  effects  of  the  small-scale  motions  are  present  in  the  above  equations  through  the  subgrid-scale 
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(SGS)  stress  tensor,  T,y,  in  the  momentum  equation  as 


T ij  =  P (m/W j  -  UiU j)  (6) 

and  the  SGS  terms  that  are  the  last  four  terms  on  the  right  hand  side  of  Eq.  (4)  as 


Qj  =  pcv(Tuj-fuj ) 

(7) 

St.-  j 

OXj 

(8) 

dui  dUj 

n  =  pd7j-pdTj 

(9) 

dui  ,  dUj 
e  =  0i^j-^d7j 

(10) 

where  Qj  is  the  SGS  heat  flux  vector.  The  terms  y,  n  and  8  were  neglected  since  only  low  Mach 
flows  were  considered  [17]. 

The  filtered  dimensionless  viscous  stress  and  heat  flux  vectors  were  approximated  by  assuming 
that  the  correlations  between  the  fluid  properties  and  the  derivatives  of  the  velocity  or  temperature 
were  weak  [22].  The  approximations  are 


G  ij  ~  Gy  — 


Re  ref 


Sij  k^^kk^i 


(ID 


and 


where  the  strain  rate  tensor  is 


-  ~  _  cp/j  dP 

RerefPr  dxj 


c  1  fdut  ,  duj 
*ij  _  2  \dxj  +  dxt 


(12) 


(13) 


To  close  the  equations,  the  SGS  stress  tensor  and  heat  flux  vector  in  the  Favre-filtered  equations 
need  be  modeled.  A  dynamic  model  proposed  for  compressible  turbulence  by  Moin  et  al.  [7]  and 
recommended  by  Lilly  [23]  was  implemented. 
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The  anisotropic  part  of  the  SGS  stress  tensor  based  on  the  grid  filter  was  modeled  as 


% ij  ~  ^Ikk&ij  —  ~ 2CrfpA2|5|  ^§ij  ^Skk&ij^J 
And  the  isotropic  part  of  SGS  stress  tensor,  xkk,  becomes 

xkk  =  2QpA2|S|2 


(14) 


(15) 


The  isotropic  part  of  SGS  stress  tensor  was  neglected  because  it  has  a  lower  order  of  magnitude  than 
the  thermodynamic  pressure  [7].  The  coefficient,  Q,  was  computed  dynamically  by 


W 

(MuMu) 


(16) 


where 


Lij  =  ptijUj 


pUiPUj 


Mu  =  -2pA2|l|(l,v  -  ^5kk&ij)  +  2A2p|.S|  (Sjj  - 


(17) 

(18) 


A  procedure  similar  to  the  modeling  of  the  SGS  stress  tensor  was  followed  to  represent  the  SGS 
heat  flux  vector.  Considering  the  modeling  for  the  eddy  diffusivity  SGS  model,  the  subgrid-scale 
heat  flux  vector  can  be  modeled  as 


n.  - 

Pr,  dxj 

cpQpA2lS|  3 f 
Pr,  dxj 

where  Pr,  is  the  turbulent  Prandtl  number  and  calculated  dynamically  as 


Pr,  =  -cpCd 


mi 

(HkFk) 


where 


(19) 


(20) 


(21) 
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Fi  =  pi’-|l|^-A2(p|S||) 


(22) 
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3  Channel  Flows  with  Rotation 


3.1  Introduction 

Turbulent  rotating  flows  are  of  considerable  interest  in  a  variety  of  industrial,  geophysical  and 
astrophysical  applications.  Examples  are  natural  flows  like  ocean  currents,  estuaries  and  atmo¬ 
spheric  boundary  layers,  and  engineering  flows  in  rotating  devices  such  as  turbines,  pumps,  com¬ 
pressors  and  cyclone  separators.  It  is  well  known  that  system  rotation  affects  both  the  mean  motion 
and  the  turbulence  statistics.  Many  rotational-induced  flow  phenomena  have  been  reviewed  by  Trit- 
ton  [4]  and  Hopfinger  [5]. 

System  rotation  influences  turbulence  in  several  ways:  for  example,  it  may  decrease  energy 
transfer  from  large  to  small  scales  or  reduce  turbulence  dissipation  and  the  decay  rate  of  turbulence 
energy.  In  rotating  channel  flow,  system  rotation  can  both  stabilize  and  destabilize  the  flow.  On 
the  unstable  side,  Coriolis  forces  enhance  turbulence  production  and  increase  the  intensity  of  tur¬ 
bulence,  while  on  the  stable  side,  Coriolis  forces  reduce  turbulence  production  and  decrease  the 
intensity  of  turbulence  (Piomelli  and  Liu  [6]). 

Most  large  eddy  simulations  (LES)  of  flows  with  heat  transfer  reported  to  date  employed  an 
incompressible  formulation  and  treated  temperature  as  a  passive  scalar.  However,  compressible 
formulations  have  been  employed  in  a  few  recent  works  (Moin  et  al.  [7]  and  Dailey  et  al.  [2])  where 
the  coupling  between  velocity  and  temperature  fields  were  considered. 

Rotating  channel  flows  have  been  investigated  experimentally  (Johnston  [8];  Han  et  al.  [9])  and 
numerically  (Piomelli  [6];  Kristofferson  and  Anderson  [10]).  However,  very  little  information  is 
available  for  rotating  channel  flow  with  heat  transfer,  especially  from  the  LES  or  DNS  community. 

The  flow  over  two  and/or  three-dimensional  obstacles  of  different  shapes  and  sizes  with  and 
without  rotation  have  been  studied  extensively  by  numerous  investigators  due  to  its  importance 
to  engineering  applications.  Among  these  are  flows  in  turbines,  pumps,  diffusers,  and  electronic 
components  (Matsubara  and  Alfreson  [11]).  In  many  of  these  applications,  enhanced  surfaces  and 
rotation  significantly  alter  the  structure  of  the  turbulence.  Han  et  al.  [12]  conducted  an  experimental 
study  to  investigate  the  effect  of  rib  geometry  on  friction  factor  and  Stanton  number  for  turbulent 
flow.  It  was  found  that  the  shape  of  the  rib  affected  the  friction  factor,  while  a  modest  effect  was 
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observed  on  the  heat  transfer. 

Bergeles  and  Athanassiadis  [13]  studied  the  influence  of  the  stream  wise  length  of  a  rib  on  reat¬ 
tachment  length  and  showed  that  a  sudden  decrease  in  reattachment  length  from  1 1  to  3  rib  heights 
was  observed  when  the  length  to  height  ratio  of  a  rib  was  greater  than  4.  Sparrow  and  Tao  [14] 
used  the  naphthalene  sublimation  technique  in  flat  rectangular  channels  of  large  aspect  ratios  with 
obstacles  situated  on  one  of  the  walls  of  the  channel  and  oriented  transversely  to  the  flow  direc¬ 
tion.  The  results  showed  a  substantial  enhancement  of  Sherwood  numbers  (Sh)  compared  with  the 
smooth- wall  duct.  Drain  and  Martin  [15]  performed  laser  doppler  velocimetry(LDV)  measurements 
of  the  fully  developed  water  flow  in  a  rectangular  duct  with  one  surface  roughened  with  a  periodic 
array  of  elements.  They  found  that  the  k- s  turbulence  model  tended  to  seriously  underestimate  the 
reattachment  length,  which  is  an  important  indicator  of  turbulence  structure. 

According  to  Wagner  et  al.  [16],  approximately  75  %  of  the  estimated  uncertainty  in  calculating 
the  heat  transfer  coefficient  was  due  to  the  temperature  measurement  error.  Furthermore,  it  can  be 
very  difficult  and  expensive  to  obtain  detailed  information  about  the  flow  distribution  in  a  ribbed 
rotating  passage  experimentally.  Large  eddy  simulation  presents  an  attractive  alternative  to  experi¬ 
ments  for  studying  details  of  such  flows.  The  goal  of  the  present  study  was  to  perform  a  large  eddy 
simulation  of  rotating  turbulent  flow  in  a  plane  channel  with  and  without  transverse  square  ribs  on 
the  walls.  Periodic  and  step  periodic  (Dailey  et  al.  [2])  boundary  conditions  were  used  at  the  inflow 
and  outflow  boundaries  since  fully  developed  channel  flows  were  considered.  This  assumption  al¬ 
lows  the  computed  domain  to  be  limited  to  the  region  between  two  adjacent  ribs  so  that  a  reasonable 
computational  grid  can  be  used. 

3.2  Numerical  Method 

In  the  present  research,  a  compressible  LES  formulation  was  used,  which  also  permited  the 
subgrid-scale  turbulent  Prandtl  number  to  be  computed  dynamically.  A  finite  volume  method  was 
used  to  numerically  solve  the  governing  equations.  The  code  used  Cartesian  hexahedral  control 
volumes,  and  solved  for  the  primitive  variables  (  p,u,v,w,T)  that  were  stored  at  the  cell  centers. 
Besides  the  physical  time  integration,  pseudo  time  iterations  were  performed  to  resolve  nonlineari¬ 
ties  in  the  algebraic  formulation  and  to  implement  the  low  Mach  number  preconditioning  using  an 
implicit  LU-SGS  scheme.  The  low  Mach  number  preconditioning  was  used  to  enable  the  compress- 
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ible  code  to  work  efficiently  at  nearly  incompressible  speeds.  The  solver  is  second-order  accurate 
in  both  space  and  time  (Dailey  and  Pletcher  [18]).  The  multiblock  code  was  parallelized  using  the 
message  passing  interface  (MPI).  The  computations  were  carried  out  on  an  IBM  SP2  (Minnesota 
Supercomputing  Institute)  using  17  processors. 

When  property  variations  are  taken  into  account,  flows  with  heating  or  cooling  do  not  attain 
a  fully-developed  state.  However,  experiments  show  that  far  downstream  of  the  entry  region,  a 
slowly  evolving  quasi-developed  state  exists.  Dailey  and  Pletcher  [18]  showed  that  a  short  section 
of  the  downstream  region  could  be  computed  in  a  “stepwise  periodic”  manner  with  the  following 
streamwise  boundary  conditions. 

pu(0,y)  =  P  u(Lx,y) 
v(0,y)  =  v(Lx,y ) 
w(0,y)  =  w(Lx,y) 

PP(0,y)  =  pP(Lx,y) 

T(0,y)  =  T(Lx,y)-ATx  (23) 

where  Lx  is  the  length  of  the  computation  domain  in  the  streamwise  direction  and  pp  is  the 
periodic  component  of  the  pressure  p  (x,  y,  z,  t)  =  fix + pp(x,y,z,t),  where  (3  is  the  average  streamwise 
pressure  gradient.  The  temperature  difference  ATX  is  computed  from  an  energy  balance  utilizing  the 
uniform  heat  flux  imposed  and  the  mass  flux.  All  primitive  variables  (p,Uj,T)  were  assumed  to  be 
periodic  in  the  span  wise  (z)  direction. 

The  conventional  Reynolds  (or  ensemble)  average  of  a  quantity  is  denoted  as  <>,  and  the  Favre 
ensemble  average  as  {},  where 

{/}  =<  P/  >  /  <  P  >  •  <24) 

A  single  prime ',  and  a  double  prime  ",  denote  the  turbulent  fluctuations  with  respect  to  the  Reynolds 
or  Favre  ensemble  average,  respectively. 

The  velocity  fluctuations  were  obtained  at  each  time  step  as 

u'i(x,y)  =  Ui(x,y)-  <  ut  >z  (x,y)  (25) 


10 


where  <>z  denotes  an  average  in  the  z  direction  only.  The  ensemble  averaged  root-mean-square 
(i rms )  values  were  subsequently  obtained  as 


'<  u'(x,y)2  > 


N, 


stat 


(26) 


where  <>  denotes  an  average  in  z  and  in  time,  and  Nstat  is  the  number  of  time  steps  used  to  com¬ 
pute  the  statistics.  The  rotating  channel  flow  simulation  began  with  a  fully-developed  channel  flow 
solution  obtained  from  previous  research,  and  needed  about  10,000  physical  time  steps  to  become 
statistically  stationary  under  the  influence  of  the  system  rotation.  The  turbulent  statistics  (e.g  urms) 
were  collected  over  the  following  10,000  physical  time  steps. 

Shear  stress  and  heat  flux  distributions  were  also  computed  as  part  of  the  turbulence  statistics. 
The  computed  shear  stress  contributions  were 


^res 

^vis 


'Zsgs 


-  <  pwV'  > 
u  du 

~  <  TTT-  > 
Re  dy 

du 

~<*Ty> 


(27) 

(28) 

(29) 


where  zres  is  the  resolvable  Reynolds  shear  stress,  xvis  is  the  viscous  shear  stress,  and  xsgs  is  the 
modeled  SGS  stress.  Similarly,  the  computed  heat  flux  contributions  were 


Qres  — 

-  <  p v"f "  > 

(30) 

Qcon  — 

.  Rcp  df  . 

RePr  dy 

(31) 

£ 

II 

/ itcv  d T 

Pr,  dy 

(32) 

where  qres  is  the  resolvable  turbulent  heat  flux,  qcon  is  the  heat  conduction,  and  qsgs  is  the  modeled 
SGS  heat  flux. 
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3.3  Results  and  Discussion 


3.3.1  Channel  Flows  without  Ribs 

The  simulations  were  run  with  a  domain  size  of  2n  x  2  x  n  with  a  grid  that  had  48  x  64  x  48 
control  volumes  in  the  x,  y  and  z  directions,  respectively.  A  grid  study  (Dailey  [2])  has  shown  that 
for  low  heating  simulations  without  rotation,  this  grid  size  provided  accurate  results  compared  to 
DNS  and  experimental  data.  The  near  wall  region  is  well-resolved.  For  the  high  heating  case,  the 
y+  (based  on  averaged  friction  velocity)  at  the  first  grid  point  near  the  wall  is  less  than  0.5. 

No-slip  velocity  and  zero  normal  pressure  gradient  boundary  conditions  were  enforced  at  the 
upper  and  lower  walls.  The  isoflux  thermal  wall  boundary  conditions  were  used  for  both  walls.  The 
dimensionless  wall  heat  flux  (qw  =  was  kept  at  2  x  10-4  for  the  low  heating  case,  and 

kept  at  2  x  10-3  and  -2  x  10-3  for  high  heating  and  high  cooling  cases,  respectively. 

Calculations  were  performed  for  four  different  cases:  no  heat  transfer,  low  heating,  high  heating, 
and  high  cooling.  The  Reynolds  number  Re  =  U{ 28)/v  (based  on  the  channel  width,  25,  and 
reference  velocity,  U)  was  5600,  and  the  rotation  number  Rob  =  |£2|(28)/t/  was  0.144.  (Note  that 
Rossby  number  is  defined  as  Ro  =  U/{2\Q\DH),  which  is  a  quarter  of  the  inverse  of  the  rotation 
number  Rob.) 

The  results  from  first  the  case  were  compared  with  the  DNS  results  (grid  96  x  97  x  128)  of 
Piomelli  and  Liu  [6]  for  the  isothermal  rotating  channel  flow  with  Reb  =  £4(28)/v  =  5700  and 
Rob  =  0.144.  Since  it  was  found  from  previous  LES  simulations  that  the  contribution  of  SGS  model 
to  turbulent  flux  is  small,  it  is  interesting  to  do  a  simulation  (coarse  grid  DNS)  without  any  model 
using  the  same  grid  resolution  as  that  of  LES  (48  x  64  x  48). 

Figure  1  shows  the  rms  distributions  (normalized  by  the  average  value  of  ux  on  two  walls) 
for  velocity  fluctuations  in  x,  y,  and  z  directions,  respectively.  The  rms  fluctuations  are  enhanced 
near  the  unstable  side  (trailing,  y  =  1.0)  but  reduced  near  the  stable  side  of  the  channel  (leading, 
y  =  -1.0).  It  is  noticeable  that  the  profile  of  vrms  becomes  a  one-peak  distribution  instead  of  the 
two-peak  distribution  obtained  without  rotation.  This  is  also  due  to  the  fact  that  the  Coriolis  force 
acts  in  the  y  direction.  These  results  agree  well  with  those  of  Piomelli[6]’s,  validating  the  current 
LES  formulation.  The  coarse  grid  DNS  data,  in  general,  slightly  overestimate  velocity  fluctuations 
compared  with  the  LES  data.  There  is  an  overshoot  for  the  coarse  grid  DNS  in  the  prediction  of 
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velocity  fluctuations  near  the  upper  (top)  wall.  The  reason  might  be  that  the  current  grid  resolution 
(48  x  64  x  48)  is  not  fine  enough  for  DNS.  LES  has  better  damping  behavior  at  the  near  wall  region 
with  the  aid  of  the  SGS  model. 

Tables  1  (a)  and  (b)  list  some  parameters  at  both  walls, 
where  the  friction  velocity  is 


Ux  — 


(33) 


A T  —  TWitop  —  TWtbol  is  the  temperature  change  across  the  normal  direction,  and  7),  is  the  bulk  tem¬ 


perature. 

T  =  S-jpuTdy 
b  fl\pudy 


Table  1  (a) 

Summary  of  the  numerical  calculations 


Case 

,tOp 

Mz,bot 

AT 

Tb 

Low  Heating 

0.0701 

0.0470 

-0.0585 

High  Heating 

0.0889 

0.0759 

-0.594 

High  Cooling 

0.0686 

0.0263 

+0.471 

(34) 


Table  1  (b) 


Case 

Nulop 

Nubot 

Cf,top 

CfM 

Low  Heating 

60.04 

20.70 

0.00551 

0.00228 

High  Heating 

66.29 

18.47 

0.00587 

0.00279 

High  Cooling 

62.14 

18.57 

0.00544 

0.00212 

The  Nusselt  numbers  and  the  friction  coefficient  are  defined  as 


Nu  = 


m 

kb 


(35) 
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(36) 


c  —  — 

f  ~  jPbUt,2 

where  h  is  the  heat  transfer  coefficient,  kb  the  bulk  thermal  conductivity,  p b  and  Ub  are  the  bulk 
density  and  velocity,  respectively. 

The  ratio  of  friction  velocity  at  the  top  wall  to  that  at  the  bottom  wall  is  approximately  2.6, 
1.5,  and  1.2  for  high  cooling,  low  heating,  and  high  heating,  respectively.  The  high  cooling  case 
reacted  most  strongly  to  the  system  rotation  as  indicated  by  the  significant  difference  between  its 
two  friction  velocities.  The  low  heating  case  has  a  relative  temperature  change  (A T / Tb )  of  5.85%, 
but  high  heating  and  high  cooling  cases  have  much  larger  temperature  change,  namely,  59.4%  and 
47.1%,  respectively.  It  is  found  that  Nusselt  number  at  the  top  wall  (or  the  unstable  side)  is  slightly 
more  than  three  times  as  large  as  that  at  the  bottom  wall  (or  the  stable  side).  Friction  coefficients  at 
the  top  wall  are  around  two  times  larger  than  those  at  the  bottom  wall. 

Figure  2  illustrates  the  streamwise  velocity  profiles  for  high  heating  (HH),  high  cooling  (HC) 
and  low  heating  (LH)  cases  in  global  and  semi-local  coordinate  systems.  The  influence  of  rotation 
is  reflected  in  the  asymmetric  distribution  of  velocities  for  all  three  cases.  The  U  velocity  profile  for 
HC  is  above  that  of  LH,  but  the  profile  for  HH  is  below  that  of  LH.  This  is,  in  part,  because  in  the 
figures  the  streamwise  velocity  has  been  normalized  by  the  average  of  uT  at  the  two  walls  and  it-  is 
significantly  influenced  by  property  variations  in  the  flow. 

Huang  et  al.  [20]  have  proposed  a  semi-local  definition  of  friction  velocity,  in  which  the  local 
density  is  used  instead  of  the  wall  density. 


When  u*  is  used  to  replace  u^,  the  resulting  semi-local  distribution  shows  that  the  three  profiles 
collapsed  toward  the  low  heating  curve. 

The  mean  temperature  profiles  are  plotted  in  global  coordinates  and  semi-local  coordinates  in 
Fig.  3.  The  Tx  and  T*  are  defined  as  below, 

Tx  =  7T^’  T?  =  (37) 

pwCpiix  p(y)CpUx 
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As  for  the  velocity  profiles,  the  high  heating  and  high  cooling  temperature  profiles  in  global  coordi¬ 
nates  departed  from  the  low  heating  results.  However,  the  spread  in  the  curves  is  greatly  decreased 
when  semi-local  coordinates  are  used.  Again,  non-symmetry  has  been  observed  in  those  averaged 
temperature  distributions. 

The  root-mean-squares  of  the  velocity  fluctuations  with  respect  to  Reynolds  ensemble  averages 
are  plotted  in  global  coordinates  in  Fig.  4.  Again,  larger  velocity  fluctuations  were  found  near 
the  unstable  side  as  a  consequence  of  rotation,  and  the  high  heating  and  cooling  values  varied 
significantly  from  the  low  heating  result.  Again,  the  values  are  normalized  by  the  average  of  uz  for 
the  two  walls. 

Figure  5  shows  the  viscous,  resolved,  and  modeled  SGS  shear  stress  distributions  normalized 
by  the  average  of  the  two  wall  shear  stresses.  The  results  indicate  that  rotation  increased  the  shear 
stress  near  the  unstable  side  but  suppressed  the  shear  stress  near  the  stable  side. 

Figure  6  shows  the  heat  conduction,  resolved  turbulent  heat  flux  and  modeled  SGS  heat  flux 
distribution  normalized  by  the  wall  heat  flux.  The  same  trends  were  observed  as  for  the  shear  stress 
distributions,  suggesting  that  heat  transfer  was  increased  near  the  unstable  side  but  decreased  near 
the  stable  side  because  of  rotation. 

3.3.2  Channel  Flows  with  Ribs 

The  geometry  and  coordinate  systems  for  the  ribbed  channel  cases  are  depicted  in  Fig.  7.  The 
ribs  were  directly  opposed  and  aligned  normal  to  the  main  flow  direction.  The  ratio  of  rib  spacing 
to  rib  height  was  10  and  rib  height  to  channel  height  was  0.1.  Calculations  were  performed  for 
three  different  cases  with  a  (72  x62  x72)  grid  size: no  rotation,  medium  rotation,  and  high  rotation. 
Uniform  heat  flux  was  applied  to  the  channel  to  investigate  the  effects  of  property  variations  with 
temperature.  The  Reynolds  number  was  5600.  The  detailed  information  is  given  in  Table  2. 


17 


Figure  7:  Schematic  of  the  computational  domain  for  the  ribbed  channel 


Table  2 

Parameters  for  three  rotational  cases 


CASE 

Rob 

I 

0.00  (No  Rotation) 

2  x  10~3 

II 

0.15  (Medium  Rotation) 

2  x  10-3 

III 

0.30  (High  Rotation) 

2  x  10“3 

Figure  8  shows  the  corresponding  rms  distributions  for  turbulent  intensities  in  the  x,  y,  and  z 
directions  at  section  A.  As  we  observed  earlier,  the  Coriolis  force  tends  to  reduce  and  enhance  the 
turbulent  intensities  near  the  stable  side,  and  unstable  side,  respectively. 

The  friction  coefficient  (C/)  and  the  Nusselt  number  ( Nu )  along  the  wall  at  sections  B  and  C  are 
shown  in  Fig.  9.  Since  the  profile  is  symmetric,  only  one  distribution  of  C/  is  shown  for  Case  I.  The 
reattachment  length  varied  with  different  rotation  numbers.  Increasing  rotation  led  to  a  decrease  (or 
increase)  in  the  reattachment  length  on  the  unstable  (stable)  side  because  the  Coriolis  force  increases 
(decreases)  turbulent  mixing  near  the  unstable  (stable)  side.  The  enhanced  mixing  on  the  unstable 
side  resulted  in  higher  levels  of  Nusselt  number  as  can  be  seen  in  Fig.  9. 
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Overall,  the  Nusselt  number  increases  rapidly  and  reaches  a  local  maximum  near  the  upstream 
comer  of  the  rib  where  high  levels  of  turbulent  intensities  are  observed.  Then,  the  Nusselt  number 
decreases  near  the  top  of  the  rib.  The  Nusselt  number  keeps  decreasing  along  the  rear  face  of  the 
rib  due  to  the  recirculation  and  increases  to  the  second  local  maximum  near  reattachment.  The  heat 
transfer  on  the  top  of  the  rib  increases  about  35  %  compared  to  stationary  (non-rotating)  case  near 
the  unstable  side,  but  decreases  about  20  %  near  the  stable  side.  Similar  order  of  magnitude  dif¬ 
ferences  compared  to  the  stationary  case  are  noticed  in  the  downstream  region.  The  flow  pattern  of 
streamlines  for  Case  I  (no  rotation)  and  Case  III  (high  rotation)  are  shown  in  Fig.  10.  The  separation 
bubbles  on  the  front  and  rear  faces  of  the  ribs  are  clearly  visible.  Due  to  the  large  increasing  pressure 
gradient  on  the  top  of  the  ribs,  a  small  separation  is  induced  (see  Fig.  9).  Behind  the  ribs,  a  rapid 
rise  of  the  static  pressure  causes  a  larger  separation  region.  For  Case  I,  the  reattachment  length  of 
the  separation  region  is  about  5h  (5  x  rib  height)  for  both  walls,  but  it  is  about  3h  near  the  unstable 
side  and  5.5h  near  the  stable  side  for  Case  III. 

3.4  Concluding  Comments-Rotating  Channels 

The  LES  results  (mean  velocity  and  velocity  fluctuations)  agree  well  with  Piomelli  [6]’s  DNS 
data  for  the  incompressible  isothermal  rotating  channel  flow  at  a  rotation  number  of  0.144. 

Turbulent  channel  flow  with  low  heating  was  calculated.  Additionally,  turbulent  channel  flows 
with  constant  wall  heating  or  cooling  rates  of  magnitudes  large  enough  to  cause  significant  variation 
in  the  temperature-dependent  fluid  properties  were  simulated.  All  those  simulations  were  performed 
under  the  influence  of  spanwise  system  rotation. 

In  general,  the  system  rotation  was  found  to  suppress  turbulent  velocity  fluctuations  and  shear 
stress  near  the  stable  side  of  the  channel,  but  enhance  those  fluctuations  and  shear  stress  near  the 
unstable  side.  Accordingly,  turbulent  temperature  fluctuations  and  turbulent  heat  flux  are  decreased 
near  the  stable  side  of  the  channel,  but  increased  near  the  unstable  side  of  the  channel.  The  stream- 
wise  and  spanwise  rms  velocity  fluctuations  can  be  seen  to  differ  by  about  a  factor  of  two  near  the 
two  walls.  The  ratio  of  wall  shear  stress  and  heat  flux  on  the  two  sides  ranged  between  a  factor  of 
2  and  3.  The  ratio  of  Nusselt  numbers  at  the  two  walls  was  very  nearly  a  factor  of  3  as  can  be  seen 
from  Table  1(b). 

The  coarse  grid  DNS  gave  nearly  as  good  results  as  LES,  except  that  it  overestimated  the  veloc- 
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Friction  Factor 


Nusselt  Number 


Figure  9:  C7  and  Nu  profiles  at  section  B  and  C  (see  Fig.  7  for  section  B  and  C) 


Case  I,  No  Rotation 


Case  III,  High  Rotation 


Figure  10:  Streamlines  for  Cases  I  and  III 
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ity  fluctuations  near  the  wall  due  to  its  lack  of  a  dissipation  mechanism.  In  addition,  using  semi-local 
coordinates  instead  of  wall  coordinates  tends  to  collapse  curves  of  high  heating  and  low  heating  into 
that  of  low  heating.  It  is  interesting,  too,  that  the  distribution  of  rms  of  V  velocity  is  changed  sig¬ 
nificantly  (from  a  two-peak  profile  to  a  one-peak  profile)  due  to  the  Coriolis  force  acting  in  the  y 
direction. 

An  investigation  was  also  conducted  on  the  effects  of  rotation  on  the  heat  transfer  in  a  ribbed 
channel  flow  using  large  eddy  simulation.  Very  consistent  results  were  obtained  indicating  that 
rotation  stabilizes  the  leading  side  and  unstabilizes  the  trailing  side.  Urms,  Vrms,  Wrms,  and  friction 
factor,  all  consistently  show  that  the  turbulence  structure  is  greatly  altered  by  the  rotation.  The 
turbulence  intensity  is  significantly  higher  near  the  unstable  side.  The  reattachment  length  also 
shows  sensitivity  to  rotation.  The  reattachment  length  on  the  stable  side  is  increased  due  to  the 
laminarization,  while  it  is  decreased  on  the  unstable  side.  Rotation  also  influences  the  heat  transfer 
in  the  flow.  With  rotation,  the  heat  transfer  was  greatly  enhanced  on  the  unstable  side  and  reduced 
on  the  stable  side.  It  is  found  that  the  temperature  profile  is  no  longer  symmetric  due  to  the  Coriolis 
force. 
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4  Large  Eddy  Simulation  of  Turbulent  Heat  Transfer  in  a  Square 
Duct  with  and  without  Rotation 


4.1  Introduction 

Turbulent  flow  inside  a  rotating  square  duct  has  many  engineering  applications  including  turbine 
blade  cooling.  In  turbulent  duct  flows,  Prandtl’s  secondary  flow  of  the  second  kind  has  a  significant 
effect  on  the  transport  of  heat  and  momentum.  This  secondary  flow  causes  distortion  of  the  isolines 
of  mean  velocity,  temperature,  and  their  fluctuations.  Heating  changes  the  structures  of  the  near  wall 
turbulence  which  in  turn  changes  the  mean  and  fluctuation  values  as  well  as  Reynolds  stress  compo¬ 
nents.  Different  heating  arrangements  have  different  effects.  With  rotation,  Coriolis  and  centrifugal 
buoyancy  forces  cause  a  more  complicated  secondary  flow  pattern  and,  as  a  consequence,  modify 
heat  transfer  coefficients.  The  objective  of  this  portion  of  the  study  is  to  expand  the  capability  of 
the  current  LES  code  to  correctly  predict  complex  turbulent  flow  phenomena  and  to  gain  a  better 
understanding  of  the  physics  of  turbulent  flow  in  rotating  passages.  In  the  following  sections,  we 
will  present  the  numerical  approach  first  and  then  the  results  of  various  cases. 

4.2  Approach 

The  numerical  formulation  used  in  the  research  is  based  on  a  finite  volume  discretization  of  the 
Favre-filtered  compressible  Navier-Stokes  equations.  A  finite  volume  LU  decomposition  scheme 
coupled  with  time  derivative  preconditioning  is  used  to  simulate  compressible  three-dimensional 
turbulent  flows  at  low  Mach  numbers.  The  whole  calculation  domain  is  divided  into  two  parts  (see 
Fig.  11).  Spatially  periodic  boundary  conditions  are  applied  to  the  first  (or,  the  inflow  generator) 
part  of  the  domain,  and  the  fully  developed  exit  velocity  profile  of  this  part  is  used  as  the  entrance 
profile  of  the  second  (or,  the  test  section)  part.  The  pressure  at  the  same  entrance  is  interpolated 
from  the  interior  domain,  however.  To  ensure  the  mass  flow  rates  of  the  two  parts  the  same,  the 
temperature  at  the  entrance  of  the  test  section  is  recalculated  according  to  the  interpolated  pressure 
so  that  the  density  is  unchanged  across  the  entrance  section.  Characteristic  out-flow  conditions  are 
applied  to  the  second  part.  This  allows  the  flow  to  develop  further  as  it  responds  to  the  heating 
condition. 
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4.2.1  LU  Decomposition  Method 


The  filtered  nondimensional  compressible  Navier-Stokes  equations  can  be  recast  in  vector  form 


au  dFi_, 

dt  a m  ’ 


where 


U=  (p,pn,pv,pw,pe)7 


The  resolved  fluxes  F;  are 


F  t  = 


P  Mf 

pWjMl  -  6,1  +  Tji 
pufi2  -  &a  +  in 
pn;n3  -  ct,3  +  X/3 

\  pen,- w/6i7  +  4  +  9»i  / 


where 


6(7  —  ■p5j7  +  ^~  (4;  ^Skk$ij)i 

Xjj  =  p(UjUj  —  Ufij ), 


e  =  cvr  +  -tiiUj, 


9(  = 


cPV 


dt 


Re,Pr  dxj 


(38) 


(39) 


(40) 


(41) 

(42) 

(43) 

(44) 
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and 


(45) 


qti  =  pcv(Tui-Tuj). 


The  source  term  S  includes  possible  body  forces.  In  the  above  equations,  the  subgrid-scale  stress 
and  heat  flux  tensors  Ty,  qtl  are  modeled  with  the  localized  dynamic  SGS  model  of  Piomelli  and 
Liu  [6].  Coupled  with  the  equation  of  state  p  =  p RT,  the  above  system  can  be  solved  by  a  LU 
decomposition  method.  For  the  sake  of  notational  simplicity,  in  the  following  text  the  overbar  for 
the  grid  scale  variables  will  be  dropped.  Equation  (38)  can  be  recast  in  terms  of  primitive  variables 
as 


[T] 


dW  3F; 
dt  +  dxi 


=  s, 


(46) 


where 

W  =  (p,u,v,w,T)T , 


(47) 


and  [T]  is  Jacobian  matrix: 


[T}  = 


au 

aw' 


Eq.  (46)  can  be  integrated  throughout  the  control  volume  Q  and  becomes 


[  [T]^-dQ  +  [  ^-dQ  =  [  SdQ. 

JQ  Ot  JQ  OXi  Jo. 


Using  the  Gauss  divergence  theorem,  the  above  equation  becomes 


[  [T]^-da+i  F i?rdS=  [  SdQ, 

Ja  dt  Jda  Ja 


(48) 


(49) 


(50) 


where  3Q  is  the  surface  surrounding  the  control  volume  Q  and  e,  is  the  unit  vector  in  x,  direction. 
When  Cartesian  hexahedral  control  volumes  are  used,  the  above  equation  can  be  approximated  in 
every  control  volume  as 

aw  6 

[T]^+^ijnjiSj  =  SQ,  (51) 

where  summation  convention  does  not  hold  for  index  j.  Fy  is  the  value  of  F,-  evaluated  at  surface 
j;  riji  is  the  projection  of  the  unit  outward  normal  Tij  of  surface  Sj  in  the  x,  direction.  Time  deriva¬ 
tive  preconditioning  developed  by  Pletcher  and  Chen  [3]  is  adopted  to  overcome  the  singularity 
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caused  by  low  Mach  number.  This  preconditioning  introduces  a  pseudo  temporal  derivative  into  the 


equation: 


[r]^-Q  +  [r]^-Q+  'ZFijnjiSj  =  SQ, 


(52) 


where  [r]  is  the  preconditioning  matrix  and  x  is  the  pseudo  temporal  variable.  This  pseudo  temporal 
derivative  term  is  discretized  as 


3W  AW 
1k~  Ax  ’ 


(53) 


where 


AW  =  W("'+1)  —  W(m) . 


(54) 


and  m  denotes  the  pseudo  time  step.  The  physical  temporal  derivative  is  evaluated  with  a  second 
order,  3 -point  finite  difference  scheme 


aw  3W("'+1)  -  2WW  +  w(',_1) 

dt  2A t 


(55) 


or  in  the  “A  fomT’as 

aw  3AW  3W(m)  -  2WW  +  w(n_1) 
~dt  2Ar +  2  At 


(56) 


where  n  denotes  the  physical  time  step.  The  flux  term  F,-  is  decomposed  into  two  parts:  the  inviscid 
part  and  another  part  F,-,  which  takes  care  of  the  viscous  and  subgrid-scale  terms.  J*}  is  treated 
implicitly  while  F,  is  treated  explicitly.  at  pseudo  time  step  m  + 1  is  approximated  as 


jj-(»'+i)  ^  j?(»») 


+  [A];  AW, 


(57) 


where 


m  ■  = 

1  jl  aw  * 


(58) 


Thus 


+  ([A],-AW)7  +  F^  =  F \J>  +  ([A],  AW);, 


("')  _  rW  . 


(59) 
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in  which  is  evaluated  with  at  face  j  (which  can  be  obtained  by  averaging  W(w)  of  two 
neighboring  control  volumes);  while  ([A], AW);  can  be  split  into 

([A], AW);  «  ([A]+AW)y_i  +  ([A]rAW)y+i,  (60) 

where  7-1  and  7+1  denote  two  neighboring  control  volumes  separated  by  face  j  and 

[A}?  =  ^[A]i±max\\[A]i\[I],  (61) 

in  which  are  the  eigenvalues  of  [A],-.  Then  every  term  in  Eq.  (52)  is  expressed  with  W(m)  and 
AW  and  the  equation  can  be  rewritten  as 

([L]  +  [D]  +  [U])AVi  =  -3.  (62) 

Then  the  preconditioned  time  accurate  governing  equation  is  solved  by  the  lower-upper  symmetric 
Gauss-Seidel  (LU-SGS)  scheme  [24], 


4.2.2  Outflow  Boundary  Treatment 


At  the  outlet  of  the  test  section,  Navier-Stokes  characteristic  boundary  conditions  are  applied. 
This  method,  which  intends  to  provide  time-accurate  boundary  conditions,  was  proposed  by  Thomp¬ 
son  [25,  26]  and  then  was  further  developed  by  Poinsot  &  Lele  [27]. 

Consider  the  characteristic  form  of  Eq.  (46)  in  x\  direction  (since  the  normal  of  the  outlet  is  in 
x\  direction  in  present  work): 


i^+^+c=#’ 


(63) 


or 

3W  d\Y 

ir'i r+W'&+c=#’  (64) 

where  the  source  term  is  neglected  and  C  are  flux  derivative  terms  excluding  Let  [sf]  = 
[T}~1  [A] i ,  then  the  above  equation  becomes 


aw  .  ^aw  .  . A 
_  +  K]_  +  [n  1c  =  o. 


(65) 
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Let 


[A]=[srvp, 


(66) 


where  [A]  is  the  diagonal  matrix  whose  elements  are  eigenvalues  of  \sf\,  and  the  rows  of  [5]  1  are 
left  eigenvectors.  For  an  ideal  gas, 

(  A-  r  \ 


[A]  = 


U\  -f-  c 
u\  —  c 

U\ 

U\ 

V  Ml  / 


(67) 


where  c  is  the  local  sound  speed  and 


(68) 


where  y  is  the  specific  heat  ratio.  Then,  after  multiplying  [5]  ',  Eq.  (65)  becomes 


ir 1  ^ +iA]ir +[s]-'pr1c = ». 


m 


8W 


(69) 


-5? 


By  employing  the  characteristic  form,  waves  with  different  velocities  can  be  determined  separately. 
At  the  boundary,  waves  leaving  the  domain  are  calculated  using  interior  points  and  one-sided  differ¬ 
ences.  Waves  propagating  into  the  domain,  however,  should  be  estimated  by  available  information 
outside  the  domain  and  also  by  examination  of  the  above  equation.  In  the  current  situation,  j£?2  is 
the  only  incoming  wave  since  its  speed  is  negative  ( u\  —  c).  By  multiplying  [5],  Eq.  (69)  becomes 


aw 

dt 


+  [SjJf  +  [T]~'C  =  0 


(70) 


A  constant  pressure  at  infinity  is  used  as  the  outlet  boundary  condition  in  the  present  work.  This  is  a 
partially-reflecting  condition.  From  Eq.  (70),  by  neglecting  the  transverse  and  viscous  terms,  which 
is  also  called  the  local  one-dimensional  inviscid  (LODI)  assumption,  we  have 


(71) 
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Let 


1  =  0,  (72) 

dt 

then 

JS?2  =  -^1.  (73> 

Since  =  [n]^,  the  matrix  [n]  is  changed  into  [n]'  according  to  Eq.  (73).  To  take  account  of  the 
effect  of  pressure  at  infinity,  Eq.  (73)  is  modified  [27]  as 


2z?2  =  -S£\  +K(p-  poo), 


(74) 


where  K  is  a  positive  constant  and  can  be  determined  by 

<j(l  —  !M2)c 

K= - 7 - > 


(75) 


in  which  a  is  equal  to  0.27  [28],  M  is  the  maximum  Mach  number  in  the  flow,  L  is  a  characteristic 
size  of  the  domain.  The  waves  at  the  outlet  boundary  thus  becomes 


(76) 


where 


1  0  1 


K(p-poo) 


0 

0 


>■ 


And  Eq.  (69)  becomes 


(77) 


[SA  [U]'^-  +  b  +  [S]-1^]"^  =  0. 

Ot  OX  \ 


(78) 
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Multiplying  by  [T][5],Eq.  (78)  becomes 


M?+M^+c'=0' 


(79) 


or 


.  ,  aw 


Sf[ 

dxi 


+c'  =  o, 


(80) 


where 


wi  =  mm’ 

C  =  C+[rp]b.  (81) 

Equation  (80)  is  of  the  same  form  of  Eq.  (63)  and  should  be  applied  at  the  outlet  surface  as  Fig.  12 
indicates.  At  the  outlet  control  volume  ni,  the  x\  direction  flux  derivative  term  in  Eq.  (52)  should 


Figure  12:  Sketch  of  the  outlet  boundary 

be  evaluated  according  to  Eq.  (63)  at  face  j  =  3  (interior  face)  and  according  to  Eq.  (80)  at  face 
j  =  1  (outlet  boundary  face),  respectively.  Note  that  &[  is  unknown  (though  we  do  know  and 
which  is  [A]j),  thus  we  evaluate  the  integral  of  the  xi  direction  flux  derivative  term  through  the 
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Table  1:  Computational  details,  cases  without  rotation 


Case 

Reb 

Grid  Size 

Tw/Tb 

1 

8100 

180x40x40 

0.0 

1.0 

2 

16200 

180x40x40 

0.0 

1.0 

3 

8100 

180x40x40 

- 

1.23 

4 

8100 

180x40x40 

2.0  x  10"3 

- 

5 

8100 

180x40x40 

- 

2.0 

6 

16200 

180x40x40 

- 

1.23 

boundary  control  volume  ni  as 

"  (r2-F<?)5,3 

+  [([A]1AW)m--([A]1AW)3]5i3 

9xi  2 

+  [([A]/IAW)1-([A],1AW)„i]513,  (82) 

where  Si  3  is  the  area  of  face  j  —  1  and  j  —  3.  The  subscript  ni  means  the  value  at  the  center  of 
control  volume  ni.  Then  the  outflow  boundary  control  volume  can  also  be  expressed  in  the  form  of 
Eq.  (62)  and  be  solved  by  the  LU-SGS  scheme. 

4.3  Results  and  Discussion 

4.3.1  Turbulent  Heat  Transfer  in  a  Square  Duct  Without  Rotation 
The  Isothermal  Duct 

To  validate  our  code,  an  isothermal  duct  case  was  simulated.  The  computational  details  are  in 
Table  1.  The  results  are  compared  with  DNS  results  [29]  and  experimental  data  [30,  31].  The  mean 
quantities  are  obtained  by  ensemble  averaging  throughout  the  test  section  in  the  x  direction  and  in 
time.  Figure  13  shows  the  streamwise  mean  velocity  along  the  wall  bisector.  The  root  mean 
square  (r.m.s.)  values  of  u,  v  and  w  along  the  wall  bisector  are  presented  in  Fig.  14.  Both  mean 
velocity  and  velocity  fluctuations  are  normalized  by  the  friction  velocity  averaged  over  four  walls. 
If  the  mean  velocity  is  normalized  by  the  local  friction  velocity,  the  profile  will  move  closer  to  the 
Nikuradse  equation  (‘log-law’).  The  friction  velocity  and  distance  in  wall  units  (dimensional)  are 
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Figure  13:  Streamwise  mean  velocity 


defined  as 


where  T„.  is  the  wall  shear  stress  and  5V  is  the  distance  to  the  wall.  Figure  15  shows  the  local  wall 
shear  stress  normalized  by  the  averaged  wall  shear  stress  along  one  wall  compared  with  the  DNS 
results  [29].  Figure  16  shows  a  comparison  between  the  measurements  of  U/Uc  and  W  / Uc  at  differ¬ 
ent  locations  [30]  with  simulation  results  where  Uc  is  the  streamwise  mean  velocity  at  the  center  of 
the  duct.  The  W  profile  indicates  the  existence  of  secondary  flow.  The  magnitude  of  the  secondary 
flow  is  quite  low  (only  about  1%  of  the  bulk  streamwise  velocity),  however,  it  has  significant  effects 
on  the  heat  and  momentum  transport.  The  turbulent  intensities  at  different  locations  are  shown  in 
Fig.  17.  The  rms  values  at  the  wall  bisector  are  very  much  like  those  of  turbulent  channel  flows 
or  boundary  layers.  However,  the  flow  behavior  near  the  comer  is  totally  different  because  of  the 
presence  of  both  walls.  This  result  is  instructive  because  we  will  compare  our  simulation  results  at 
wall  bisector  with  DNS  results  of  a  channel  later.  Very  good  agreement  with  DNS  and  experimental 
results  was  obtained. 
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The  Heated  Duct 

The  domain  for  the  heated  duct  cases  are  the  same  as  that  of  the  isothermal  cases  which  is 
shown  in  Fig.  11.  Heat  is  applied  to  all  four  walls  of  the  test  section.  Two  different  kinds  of 
heating  are  considered:  constant  wall  temperature  (Cases  3,  5  and  6)  and  constant  wall  heat  flux 
(Case  4).  The  computational  details  are  in  Table  1.  Case  5  provides  the  highest  heating  level  in 
these  cases  while  Case  4  imposes  a  higher  heating  level  than  Case  3.  The  ratio  of  average  wall 
temperature  to  the  bulk  temperature  near  the  outlet  for  Case  3  is  1.13, 1.45  for  Case  4,  1.53  for  Case 
5  and  1.14  for  Case  6.  Because  the  heated  duct  flows  are  non  homogeneous  in  the  x  direction, 
the  mean  values  of  the  statistical  quantities  are  shown  mostly  for  a  stream  wise  location  near  the 
outlet  (except  the  streamwise  Nusselt  number  distribution)  and  are  obtained  by  averaging  only  in 
time.  Figure  18  shows  the  streamwise  mean  velocity  profiles  at  the  wall  bisector.  Here  the  mean 
velocity  is  normalized  by  the  local  friction  velocity.  For  the  same  Reynolds  number,  the  higher 
heating  velocity  profile  departs  more  from  the  log-law.  The  mean  temperature  profiles  are  plotted 
in  wall  coordinates  in  Fig.  19.  The  mean  temperature  difference  is  given  as 


/ 


Figure  14:  Turbulent  intensities 


(84) 
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y/Di, 

Figure  15:  Local  wall  shear  stress 

Also  the  higher  heating  case  departs  from  the  thermal  log-law  more.  This  trend  is  consistent  with 
what  is  observed  in  channel  flow  with  heat  transfer  under  variable  property  conditions  [2].  Figure 
20  shows  the  wall  shear  stress  distributions  normalized  by  the  averaged  wall  shear  stress.  In  the 
isothermal  case,  there  is  a  local  maximum  stress  at  the  mid-wall.  Heating  tends  to  suppress  this 
mid-wall  maximum.  This  is  consistent  with  the  results  of  one-side  heated  square  duct  cases  [32]. 
Figures  21  and  22  show  the  variation  of  heat  flux  and  temperature  around  the  duct  perimeter  for 
Cases  3  and  4,  respectively.  The  local  heat  flux  is  normalized  by  the  average  wall  heat  flux  and  the 
nondimensional  temperature  is  defined  as 

8=  Tw~(Tw\  (85) 

where  (Tw)  is  the  average  wall  temperature.  The  LES  results  of  Pallares  et  al.  [33]  and  experimental 
data  of  Brundrett  et  al.  [34]  are  shown  for  comparison  purposes.  Our  results  agree  with  the  LES 
results  well  and  the  discrepancy  with  experiment  is  probably  attributable  to  the  high  Reynolds  num¬ 
ber  the  experiment  employed  (Re  =  75000)  which  tends  to  smooth  the  heat  flux  distribution.  The 
turbulent  intensities  at  the  wall  bisector  are  shown  in  Figs.  23,  24  and  25.  DNS  results  of  a  heated 
channel  [35,  36]  and  a  thermal  boundary  layer  [37]  are  also  shown  for  comparison.  The  velocity 
fluctuations  of  the  low  heating  case  agree  with  DNS  results  well.  Higher  heating  decreases  velocity 
fluctuations  due  to  the  increased  thickness  of  the  viscous  sublayer.  The  temperature  fluctuations  of 
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z/h=0. 1 


Figure  16:  Mean  velocity  at  different  locations - ,  Case  1;  o,  Experiment  (Cheesewright 

et  al.  1990). 


Case  3  are  over  predicted  by  the  present  simulations,  especially  the  peak  value.  The  reason  is  still 
not  clear.  However,  the  location  where  the  maximum  temperature  fluctuation  occurs  for  Case  5  is  in 
good  agreement  with  DNS  data,  which  is  y+  «  18.  The  temperature  fluctuations  for  the  isoflux  wall 
case  agrees  with  DNS  results  very  well.  As  the  wall  is  approached,  the  rms  temperature  fluctuation 
becomes  a  constant  2.0  when  very  a  thin  wall  is  considered  [38].  The  distribution  of  streamwise 
Nusselt  number  is  shown  in  Fig.  26.  The  Nusselt  number  is  normalized  by  the  fully  developed  (both 
velocity  and  temperature  profile)  Nusselt  number  Nu/j.  In  the  present  work,  the  empirical  equation 

Nufd  =  0.023 Re°b%Pr0AA-°-55  (  86) 


was  used  to  estimate  the  reference  Nusselt  number.  The  series  solution  [39]  and  experimental  results 
of  Sparrow  et  al.  [40]  are  also  shown  for  comparison  purposes.  Good  agreement  has  been  obtained. 
Finally,  the  profiles  of  mean  and  rms  velocities  for  a  heated  duct  at  different  locations  are  shown  in 
Figs.  27  and  28.  Comparing  with  the  isothermal  duct,  only  the  secondary  flow  shows  some  obvious 
differences,  especially  in  the  vicinity  of  the  duct  middle  plane.  This  is  because  of  the  enhanced 
ejection  around  the  wall  bisector  which  is  caused  by  wall  heating.  This  is  displayed  on  Fig.  29. 
This  kind  of  phenomenon  has  also  been  confirmed  by  an  experimental  study  [41]. 
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Figure  17:  Turbulent  intensities  at  different  locations  - ,  Case  1;  o,  Experiment 

(Cheesewright  et  al.  1990). 


4.3.2  Turbulent  Heat  Transfer  in  a  Square  Duct  With  Rotation 


In  a  rotating  frame,  the  flow  will  feel  both  Coriolis  and  centrifugal  buoyancy  forces.  The  source 
term  in  Eq.  38  becomes 

(  0  \ 

2pRou2  — 

S  =  |  —2pRou\  I  ,  (87) 

0 
0 


when  the  rotation  is  in  z  direction  and  x  coordinate  is  aligned  with  the  flow  direction  as  Fig.  1 1 
shows.  The  rotation  number  Ro  is  Ro  =  jj^  where  Q  is  angular  velocity  and  U[,  is  the  bulk  velocity. 
The  Grashof  number,  Gr,  is  defined  as  Gr  =  _P where  rm  is  the  mean  rotating 
radius.  Notice  that  rm  can  be  negative  if  the  fluid  flows  inward  toward  the  axis  of  rotation.  Thus,  the 
Grashof  number  is  negative  for  outward  flow  and  positive  for  inward  flow  under  heating  with  the 
current  Grashof  number  definition.  £  is  equal  to  P(TV),  —  1], ) .  The  computational  details  are  in  Table 
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Figure  20:  Local  wall  shear  stress 
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Figure  21:  Local  wall  heat  flux 
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Figure  27:  Mean  velocity  at  different  locations 


,  Case  1; 


o,  Case  3; 


A,  Case  4 


Figure  29:  Instantaneous  transversal  vector  field 


A,  Case 
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Figure  34:  Secondary  velocity 


Figure  35:  Turbulent  intensities 


Two  isothermal  rotating  duct  cases  (Cases  1  and  2)  were  used  to  validate  our  code.  The  results  of 
mean  profile  and  velocity  fluctuations  at  wall  bisector  B  (see  Fig.  11)  are  compared  with  DNS  results 
[42]  (see  Figs.  30-35).  Very  good  agreement  has  been  obtained.  The  secondary  velocity  vector  and 
streamwise  velocity  contours  are  also  shown  (see  Fig.  36  to  Fig.  38).  It  can  be  seen  from  the  above 
results  that  the  Coriolis  force  generates  persistent  secondary  flows  as  well  as  increase/decrease  in 
the  mean  shear  stress  (^)  and  turbulent  intensities  at  the  unstable/stable  side,  respectively.  The 
modification  of  the  mean  shear  stress  is,  from  another  viewpoint,  the  shift  of  the  axial  flow  toward 
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Figure  36:  Secondary  flow  pattern  and  streamwise  velocity  distribution 


the  unstable  side.  When  the  rotation  number  is  high  (see  Fig.  38),  the  profile  of  U  becomes  uniform 
in  the  z  direction,  which  is  known  as  the  Taylor-Proudman  effect.  This  phenomenon  is  a  result  of  the 
balance  between  pressure  gradient  and  Coriolis  force,  which  has  no  component  in  the  z  direction. 

Then  heated  rotating  duct  cases  were  simulated.  Two  heating  conditions,  constant  wall  temper¬ 
ature  and  constant  wall  heat  flux,  were  considered.  Due  to  rotation,  the  secondary  flow  transports 
the  central  cold  fluid  to  the  unstable  side  which  results  in  relatively  lower  local  temperatures  at  the 
unstable  side  and  higher  temperatures  at  stable  side.  Density  differences  arise  as  a  consequence  of 
these  local  temperature  differences.  The  buoyancy  contributes  in  the  aiding/opposing  directions  to 
the  mean  flow  at  the  unstable/stable  sides,  respectively,  if  the  Grashof  number  is  negative  (  outward 
flow),  which  helps  in  shifting  the  axial  flow  toward  the  unstable  side  caused  by  Coriolis  force,  and 
the  situation  is  reversed  for  inward  flow  (positive  Grashof  number).  This  effect  can  be  seen  from 
Figs.  39  and  40  which  show  mean  streamwise  velocity  profiles  at  the  wall  bisectors  of  the  forced 
and  mixed  convection  cases.  Also  we  can  see  that  the  buoyancy  force  causes  the  flow  to  separate 
at  the  stable  side  when  the  negative  Grashof  number  is  large  enough  (Case  7)  in  magnitude.  Such  a 
mechanism  is  depicted  in  Fig.  43.  The  mean  temperature  profiles  at  the  wall  bisectors  of  different 
heated  cases  are  shown  in  Figs.  41  and  42.  The  nondimensional  temperature  is  defined  as 
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Figure  37:  Secondary  flow  pattern  and  streamwise  velocity  distribution 

where  (Tw)  is  the  average  wall  temperature  and  (Tz)  is  the  average  friction  temperature.  In  com¬ 
parison  of  the  buoyancy- free  flows,  a  significant  temperature  rise  appears  near  the  stable  side  for 
outward  flows  due  to  the  slower  reversed  fluid  movement  caused  by  the  opposing  buoyancy.  Such 
a  situation  always  reduces  the  heat  transfer  on  the  stable  wall.  However,  this  deterioration  is  some¬ 
what  compensated  for  by  the  enhanced  turbulence  level.  Similar  results  have  also  been  reported  by 
Hwang  et  al.  [46].  Steeper  temperature  profiles  appear  near  the  stable  wall  for  inward  flows. 

As  mentioned  above,  an  important  effect  of  rotation  on  a  turbulent  duct  flow  is  the  stabiliza¬ 
tion/destabilization  of  turbulence  on  walls  normal  to  the  rotation  axis.  This  effect  is  a  result  of 
the  modification  of  mean  shear  stress  caused  by  the  Coriolis  force.  The  buoyancy  force  influences 
turbulent  intensity  in  two  ways:  it  modifies  the  mean  shear  further  and  takes  part  in  the  turbulent 
production  term.  The  mean  velocity  fluctuation  profiles  along  wall  bisector  B  are  shown  in  Figs. 
44  and  45.  In  the  outward  flow,  the  turbulent  intensity  near  the  stable  side  is  greatly  enhanced  by 
the  buoyancy  This  is  because  the  reversed  flow  caused  by  buoyancy  gives  rise  to  a  strong  mean 
shear  gradient  near  the  stable  wall.  At  the  unstable  side,  there  is  also  some  turbulence  intensity  in¬ 
crease  because  of  the  growing  mean  shear  gradient  caused  by  the  aiding  buoyancy.  In  inward  flow, 
the  mean  velocity  fluctuation  decreases/increases  slightly  at  the  unstable/stable  side  with  the  reduc¬ 
tion/augmentation  of  mean  shear  gradients.  The  temperature  fluctuation  profiles  along  wall  bisector 
B  of  different  cases  are  shown  in  Fig.  46  and  Fig.  47.  Rotation  decreases  the  temperature  fluctu- 
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Figure  38:  Secondary  flow  pattern  and  streamwise  velocity  distribution 

ations  near  the  stable  wall  because  of  the  suppression  of  the  turbulent  heat  fluxes.  The  buoyancy 
for  the  outward  flow  (negative  Grashof  number)  considerably  increases  the  temperature  fluctuation 
intensity  at  the  stable  wall.  The  temperature  fluctuation  near  the  stable  wall  for  inward  flows  also 
increases  somewhat  because  of  the  steeper  temperature  gradient  caused  by  the  aiding  buoyancy. 

The  heat  transfer  coefficient  is  strongly  affected  by  rotation.  The  Coriolis  force  shifts  the  axial 
flow  toward  the  unstable  side  which  produces  a  thinner/thicker  boundary  layer  at  the  unstable/stable 
wall.  Due  to  the  thickness  change  of  the  boundary  layers,  the  heat  transfer  coefficient  at  the  un¬ 
stable/stable  wall  increases/decreases.  The  buoyancy  force  aids  this  tendency  for  outward  flows 
(negative  Grashof  number)  and  hinders  such  tendency  for  inward  flows  (positive  Grashof  number). 
This  can  be  seen  from  Fig.  49.  The  comparison  between  numerical  simulation  and  experiments 
[45]  is  also  given  (Figs.  48  and  49).  The  agreement  is  fairly  good  and  the  discrepancy  may  be 
attributed  to  the  differences  between  the  numerical  and  experimental  settings.  The  heat  transfer  co¬ 
efficient  profiles  for  different  rotation  numbers  and  Reynolds  numbers  are  given  in  Fig.  50  and  51. 
With  the  same  rotation  number  and  Reynolds  number,  the  isoflux  cases  show  higher  heat  transfer 
rates  on  both  unstable  and  stable  walls.  With  the  same  heating  condition  and  Reynolds  number,  the 
differences  between  heat  transfer  coefficients  on  unstable  and  stable  walls  increase  with  increasing 
rotation  number.  With  the  same  heating  condition  and  rotation  number,  the  Reynolds  number  has 
little  impact  on  the  local  heat  transfer  coefficients,  which  agrees  well  with  the  physical  meaning  of 
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Figure  39:  Streamwise  mean  velocity 


rotation  number  (which  is  a  measurement  of  the  relative  strength  of  the  Coriolis  force  compared  to 
the  inertia  force)  and  experimental  results  [44]. 


54 


Figure  41:  Mean  temperature 


Figure  42:  Mean  temperature 


Figure  43:  Sketch  of  the  flow  separation  mechanism  at  the  leading  wall  when  Gr  <  0 


Figure  44:  Turbulent  intensity:  velocity 


Figure  45:  Turbulent  intensity:  velocity 


Figure  47:  Turbulent  intensity:  temperature 


Figure  48:  Streamwise  Nusselt  number  distribution 
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Figure  49:  Streamwise  Nusselt  number  distribution 


x/D 


Figure  50:  Effect  of  rotation  number  on  streamwise  Nusselt  number  distribution 
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Figure  5 1 :  Effect  of  Reynolds  number  on  streamwise  Nusselt  number  distribution 
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5  External  Turbulent  Boundary  Layers 


5.1  Introduction 

The  successful  simulation  of  film  cooling  of  a  turbulent  flat  plate  flow  subjected  to  free  stream 
turbulence  relies  on  the  development  of  the  following  major  techniques.  First,  a  robust  and  efficient 
numerical  solver  designed  for  compressible  flows  is  needed.  Second,  a  technique  is  needed  to 
generate  appropriate  inflow  or  starting  conditions  (including  the  freestream  turbulence)  that  results 
in  the  development  of  a  flow  that  captures  the  correct  first,  second  and  even  higher  order  statistics. 
Third,  it  must  be  possible  to  accurately  resolve  the  geometrically  complex  domain  for  the  film 
cooling  configuration.  All  of  these  subtasks  are  challenging.  In  this  project  we  attempted  to  make 
several  advances  to  this  technology  that  are  discussed  in  the  following  sections.  First  we  focus  on 
the  external  boundary  layer. 

5.2  Numerical  Scheme 

In  order  to  simulate  the  compressible  turbulent  boundary  layer,  we  formulated  a  new  fractional 
step  method.  This  method  combines  a  preconditioning  technique  with  a  factorization  treatment  ror 
the  Jacobian  matrix  so  that  the  eigenvalues  would  cluster  near  one  regardless  of  the  Mach  num¬ 
ber,  and  the  threefold  coupling  among  temperature,  pressure  and  velocities,  from  which  the  major 
numerical  difficulties  arise,  would  be  decoupled.  The  applied  factorization  enables  the  resolved 
submatrices  to  be  positive  definite  which  guarantees  the  convergence  of  the  iterative  numerical  pro¬ 
cedure.  The  numerical  formulation  used  in  the  research  is  based  on  a  finite  volume  discreatiztion 
of  the  Favre-filtered  compressible  Navier-Stokes  equations.  Unlike  the  traditional  fractional  step 
method  [47],  which  was  designed  for  incompressible  flows,  this  method  offers  a  robust  way  to 
simulate  compressible  flows. 

5.2.1  Quasi-Newton  Iteration 

The  best-known  and  widely  applied  method  for  dealing  with  nonlinear  equations  is  the  Newton 
method.  There  are  two  main  considerations  for  the  Newton  like  methods,  local  convergence  and 
global  convergence.  When  the  initial  guess  is  sufficiently  close  to  the  solution,  the  Newton  method, 
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e.g.  Kantorovich  method  [48],  has  a  quadratic  convergence  rate.  Its  convergence  strictly  depends 
on  an  appropriate  initial  guess,  but  such  a  guess  is  almost  impossible  for  the  simulation  of  complex 
flows.  Alternatively,  the  quasi-Newton  method,  which  offers  a  global  convergence  advantage,  en¬ 
ables  us  to  calculate  the  complex  flows  with  a  rough  initial  guess.  But  the  price  is  the  convergence 
rate.  Unlike  the  local  method,  the  quasi-Newton  method  has  a  superlinear  convergence  rate  [48]. 
Superlinear  means  that  the  decay  of  the  residual  error  is  no  less  than  a  power  function.  For  the  sim¬ 
ulation  of  complex  flows  involving  turbulence,  superlinear  should  be  a  satisfy  convergence  speed. 
In  the  quasi-Newton  method,  the  Jacobian  matrix  will  be  adjusted  by  a  factorization  technique.  In 
this  section,  we  present  a  new  factorization  technique  to  achieve  the  global  convergence. 

Representing  the  Newton-like  method  applied  to  the  governing  equations  by 


JSv  =  -F  (89) 

where  J  is  the  Jacobian  matrix.  Usually,  for  the  discretized  Navier-Stokes  equations,  such  a  Jacobian 
matrix  may  not  be  diagonally  dominant,  and  the  condition  number  may  be  very  large,  unless  the 
very  small  time  step  is  applied.  But  even  so,  in  incompressible  flows,  the  Jacobian  matrix  still  may 
not  be  diagonally  dominant.  Therefore,  the  Jacobian  must  be  adjusted.  The  main  idea  of  the  present 
approach  is  to  search  for  a  left  multiplier  P  and  a  right  multiplier  M  so  that  the  updated  solver 
becomes 

PJMM~x&v  =  -PF  (90) 

Actually,  the  left  multiplier  is  essentially  a  preconditioning  matrix,  and  the  right  multiplier  imple¬ 
ments  the  factorization. 
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5.2.2  Preconditioning 


The  preconditioning  matrix  proposed  by  Pletcher  and  Chen  [3]  is  utilized.  It  is 

/ 


P  = 


JL. 

RCy 

TU 

P 
£  £ 
_TV_ 

d. 

P 

fw 

P 


u 

RCV 

d 

T 

P 

0 

0 

_  fu 

PCy 


y 

RCy 

0 

i 

p 


Tv 

PCy 


w 

RCy 

0 

0 

I 

a 

p 

TW 

PCy 


1  \ 

RCy 

0 
0 
0 

k  ) 


where 


*2  i2 

(17  +  F  +W  ) 


Multiplying  both  side  of  Eqs.  (2),  (3)  and  (4)  with  this  preconditioning  matrix,  yields 
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where  b=j  +  £-,e=™,  and  vector  /  =  PF.  For  simplification,  it  can  be  represented  as 
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where  G  is  the  gradient  operator,  D  is  the  divergence  operator,  a  =  8V,  = 
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and 

%i  0  0 

c=  o  ^  0 

V  0  0  K 

Here,  ^  is  a  material  derivative  and  /  is  the  identity  matrix.  Such  a  preconditioning  matrix  offers 
two  advantages.  First,  it  transfers  the  large  off-diagonal  terms  into  diagonal  part  so  that  the  nu¬ 
merical  scheme  becomes  more  robust  regardless  of  how  small  the  Mach  number  is.  Second,  this 
preconditioning  technique  enables  us  to  decouple  the  temperature  and  velocities.  Clearly,  the  three¬ 
fold  correlations  among  temperature,  pressure  and  velocities  produce  some  numerical  difficulties 
especially  in  the  low  Mach  flows. 

According  to  a  perturbation  analysis, 

8P  =  RpST  +  RTdp 

Under  the  reference  system  applied  in  this  work,  R  =  yMnlfM^f  •  This  implies  that  a  small  perturbation 
on  temperature  will  be  greatly  magnified  by  R,  and  be  brought  into  the  momentum  equations  by  8 P. 
Moreover,  the  eigenvalues  after  preconditioning  are  U  and 

[U  ( 1  +  R)  ±  y/> U~2  (R  -  1  )2  +  4Ru2}  /  (2 R)  (94) 

where  R  —  ^  and  a  is  non-dimensional  acoustic  speed,  a  =  i/y RT.  Obviously,  those  eigenvalues 
cluster  near  1 .  Before  preconditioning,  they  are  U  and  U  ±  a.  But  it  is  still  possible  to  have  zero  or 
negative  eigenvalues  in  the  solver  regardless  of  the  Mach  number.  In  subsonic  flows,  U  <  cc,  then 
U(l  +  R)  <  \JU2 {R  —  1  )2  +  4Ra2 .  In  supersonic  flows  near  the  wall  region,  the  local  velocity  U 
must  be  less  than  the  acoustic  speed  a.  In  this  region,  U(  1  +  R)  <  y/U2(R—  l)2  +  4/?a2  also.  With 
negative  eigenvalues  the  iteratilve  system  is  not  positive  definite.  Such  non  positive  eigenvalues 
can  lead  to  difficulty  with  the  quasi-Newton  solver,  and  even  cause  the  divergence  of  the  iterative 
scheme.  Therefore,  a  proper  factorization  technique  needs  to  be  developed. 
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5.2.3  Fractional  Step  Factorization 

Denote  the  matrix  M~x  in  Eq.(90)  as 


AT1  = 
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where  /  is  an  identity  matrix.  Hence,  the  Jacobian  matrix  in  Eq.  (91)  can  be  decomposed  into 
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Essentially,  this  decomposition  is  a  fractional  step  procedure.  Its  efficiency  arises  from  decou¬ 
pling  the  velocities,  pressure  and  temperature. 

5.2.4  Numerical  Procedure 

Hence,  the  full  discrete  scheme  can  be  represented  as 
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As  mentioned  above,  such  a  decomposition,  essentially,  is  a  fractional  step  procedure.  It  pro¬ 
vides  an  efficient  way  to  decouple  velocities,  pressure  and  temperature.  The  corresponding  numer¬ 
ical  procedures  are: 

1.  Resolve 

cm  =  -Si  (98) 


for  8Wi  using  an  iterative  scheme. 
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2.  Substitute  SWj  into  equations 


Sp'  = 

—f\  —  bD'SWi 

(99) 

ST'  = 

(-f5-B'SWi)/E' 

(100) 

to  obtain  8 p'  and  ST' 

3.  Resolve  the  parabolic  equation 


(a'  -  AtbD'eG)8P  =  Sp1 


(101) 


for  SP. 

4.  Obtain  velocities  by 

Sn,  =  8W,  —  C~1eGSP  (102) 

5.  Obtain  the  temperature  correction  by 

8r  =  87’/  +  £“1P(r1eG8P  (103) 

This  is  a  5  step  procedure.  Except  steps  1  and  3,  all  of  the  other  steps  can  be  solved  non- 
iteratively.  In  step  1,  C’  is  discretization  of  C;  therefore,  a  small  enough  time  step  At  will  provide  the 
diagonal  dominance,  which  should  guarantee  the  convergence  of  step  1.  In  the  traditional  factional 
step  method,  the  equation  to  be  solved  in  step  3  is  a  Poisson  equation 

AS  P  =  RHS  (104) 


But  in  compressible  flow,  the  equation  is  not  a  Poisson  equation,  but  a  parabolic  equation.  Recall 

fit  A 

that  b  =  ^  +  £-,e  =  jr  and  a  —  Equation  (101)  is  the  discretization  of  following  equation: 
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The  above  equation  can  be  written  as 


r)RP  TR 

ASP  -  (VS P)  •  (Vy )  +  Cl  (u  ■  V)SP  =  c2SP' 


(106) 


where  ci  = 


cv _ and  co  —  rc^ 

P(R+Cv)At  2  P{R+Cv)At ' 


By  the  factorization  above,  the  parts  which  need  to  be  solved  iteratively  are  step  1  and  step  3. 
Obviously,  the  associated  Jacobian  matrices  are  positive  definite,  and  diagonally  dominant.  There¬ 
fore,  they  can  be  solved  by  iterative  scheme  efficiently. 


5.3  A  Generation  Technique  for  the  Ttorbulence  in  Turbulent  Boundary  Layers 

In  order  to  generate  turbulence,  we  developed  a  dynamic  adjustment  technique  for  the  inflow 
conditions  to  establish  a  turbulent  boundary  layer  [51].  In  1988,  Spalart  [52]  introduced  a  scaled 
periodic  boundary  condition  method  to  produce  the  inlet  profile  of  spatial  evolving  turbulent  bound¬ 
ary  layers.  Lund  et  al.  [53]  further  developed  this  concept.  In  their  implementation,  instantaneous 
profiles  at  a  specific  station  were  recycled  to  the  inlet  at  each  numerical  step  after  rescaling.  This 
rescaling  was  based  on  the  similarity  laws  of  the  boundary  layer,  the  law  of  the  wall  in  the  inner 
part  and  the  defect  law  in  the  outer  part  of  boundary  layer.  Kong  et  a/.  [37]  applied  Lund  et  al.  ’s 
idea  to  temperature,  and  formulated  an  inlet  generator  for  turbulent  thermal  boundary  layers.  Their 
simulation  results  showed  that  this  treatment  also  works  well  for  the  energy  equation  when  fluid 
properties  can  be  assumed  to  be  constant. 

For  a  simulation  not  stating  from  a  proper  initial  condition,  the  above  rescaling  treatments  may 
lead  to  a  decrease  of  skin  friction  [54]  with  increasing  time.  Since  the  mean  profiles  are  adjusted 
with  the  fluctuations  synchronously,  the  flow  keeps  transferring  its  energy  to  small  structures,  and 
most  of  this  energy  eventually  goes  to  the  very  small  eddies  and  is  dissipated.  In  the  early  part  of  the 
simulation,  if  the  large  eddy  structures  have  not  been  correctly  produced,  the  flows  will  not  contain 
enough  energy  to  sustain  the  turbulence.  Finally,  the  skin  friction  will  drop  down  to  the  value 
characteristic  of  laminar  flows.  To  avoid  this,  Lund  et  al.  [53]  suggested  making  a  correction  to  the 
solved  velocities  during  the  early  part  of  simulation.  Spille-Kohoff  and  Kaltenbach  [54]  suggested 
adding  a  source  term  to  the  resolved  equation  based  on  the  desired  Reynolds  stress. 

For  solving  this  problem,  we  suggest  replacing  the  traditional  rescaling  method  with  a  dynamic 
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Figure  52:  The  configuration  for  a  turbulent  boundary  layer  and  the  recycling  downstream  station, 
which  is  dynamically  selected  by  Eq.  (107) 


adjustment  recycling  method.  This  idea  come  from  the  observation  that  the  structures  in  a  turbulent 
boundary  layer  convect  downstream  with  a  speed  proportional  to  the  streamwise  mean  velocity.  The 
proposed  recycling  method,  namely  the  dynamic  adjustment  method, can  be  represented  by 


Xrec  =  X\  +mm(Xtag-X\,aUbmax(0,(t -to)))  (107) 

where  Xrec  is  recycling  location.  Xtag  is  the  desired  recycling  station  when  the  numerical  domain  is 
completely  filled  by  the  turbulent  structures  produced  by  inflow.condition,  X\  is  the  starting  location, 
t  is  non-dimensional  time,  Ub  is  the  bulk  velocity  and  to  is  the  time  in  which  the  leading  edge  of  the 
convected  flow  generated  at  the  inlet  can  reach  station  X\.  Equation  (107  means  that  the  recycle 
plane  stays  at  station  X\  from  t  =  0  up  to  t  —  to,  when  the  convective  structures  generated  by  the 
inlet  condition  is  supposed  to  have  passed  though  station  Xi  except  for  the  viscous  sublayer,  and 
then  the  recycle  plane  moves  downstream  with  the  speed  a  Ub  until  it  reaches  the  location  of  the 
desired  recycle  station  Xtag.  After  that,  the  recycle  plane  will  stay  fixed  for  the  the  rest  of  the 
calculation.  This  dynamic  adjustment  procedure  is  designed  to  keep  the  recycle  plane  inside  the 
region  influenced  by  the  inlet  conditions  to  achieve  steady  first  and  second  order  statistics.  The 
reason  for  this  is  that  if  the  initial  condition  is  not  proper,  i.e.  the  correct  large  eddies  do  not  appear 
in  the  flow,  the  mean  profile  and  velocity  rms  profiles  may  be  deformed,  and  similarity  law  cannot 
be  applied.  However,the  similarity  laws  are  the  basis  of  rescaling  method.  Coupling  of  initial  and 
inlet  boundary  condition  may  lead  to  the  decay  of  skin  friction.  Therefore,  Lund  et  al.  suggested  a 
correction  to  velocity  profiles  during  the  early  part  of  simulation.  This  kind  of  correction  may  create 
some  convergence  problems.  And  the  physical  basis  for  this  is  not  clear.  In  our  work,  we  propose  a 
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dynamic  adjustment  as  a  remedy  for  this  problem. 

As  part  of  the  validation  procedure  a  case  was  computed  for  a  Reynolds  number  based  on  the 
displacement  thickness,  Red,  of  2000  using  a  numerical  domain  X  x  Y  xZ  —  848,*  x  308,*  x  19.28,*, 
and  a  numerical  grid  of  nt  x  tij  x  nk  =  280  x  80  x  120,  and  a  Mach  number  of  0.06.  The  time  step 
was  0.2.  A  dynamic  SGS  model  was  applied. 


Figure  53:  Comparison  of  streamwise  mean  profile  in  a  turbulent  boundary  layer  in  the  order  of 
increasing  streamwise  direction;  the  units  of  X  is  one  initial  displacement  thickness 


Figure  54:  Comparison  of  Urms  profile  in  a  turbulent  boundary  layer  in  the  order  of  increasing 
streamwise  direction;  the  units  of  X  is  one  initial  displacement  thickness 

The  streamwise  mean  profiles  at  different  streamwise  stations  are  compared  in  Fig.  53.  And 
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the  Urms  profiles  are  shown  in  Fig.  54  in  the  order  of  increasing  streamwise  direction.  From  those 
figures,  we  observe  that,  with  this  dynamic  inlet  generation  technique,  the  inlet  buffer  zone  becomes 
fairly  short,  and  the  influence  of  buffer  zone  on  the  first  and  second  order  statistics  is  very  small. 

5.4  Simulation  of  a  Turbulent  Boundary  Layer  Subjected  to  Free  Stream  Turbu¬ 
lence 

For  the  turbulent  boundary  layers  subjected  to  free-stream  turbulence,  the  transported  energy 
from  free  stream  to  boundary  layer  will  increase  skin  friction  (Blair  [55],  Hancock  and  Brad¬ 
shaw  [56]).  This  relation  was  formulated  as 

AC/  Tue 

^a2  +  L/8 

This  formula  implies  that  a  smaller  turbulent  free  stream  length  scale  makes  a  larger  contribution  to 
skin  friction  compared  with  a  larger  turbulent  free  stream  length  scale.  Since  shear  is  not  large  in 
the  free  stream,  there  is  not  much  turbulent  production  there;  therefore,  the  main  feature  in  the  free 
stream  is  dissipation.  Length  scale  is  a  parameter  related  to  the  intensity  of  dissipation.  Usually, 
grid  turbulence  has  a  length  scale  range  from  108j  to  505,/  (Barrett  and  Hollingsworth  [57],  [58]). 
In  this  range,  the  dissipation  of  free  stream  turbulence  is  very  strong.  Therefore,  the  mean  profile 
in  the  wake  region  of  a  turbulent  boundary  does  not  change  slowly  in  the  downstream  direction  as 
is  the  case  without  free  stream  turbulence.  In  the  other  words,  the  defect  law  does  not  apply  with 
free  stream  turbulence.  Those  features  are  extremely  interesting,  but  very  difficult  to  capture  in  the 
numerical  simulations. 

In  order  to  describe  this  phenomena,  we  decompose  the  velocity  into 

UlMter(ytni,)  =  u  +  Hx,T\)up(x,r\,Z,t)  (108) 

U°nlter(ainlt )  =  (Y I  UrecyiU-mll )  +  Y2  U'recy{CLinlt  ,Z,t))  +  (1  -  Yl  )U„f  (109) 

Ufn7( 0W,)=Y3(  X  Weighty'^) +  Uref  (110) 

0  <;<m 
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where  the  inner  part  ranges  from  0  to  30  in  the  y+  sense.  Spalart’s  [52]  inlet  solver  was  applied 
for  the  inner  part  because  the  correlation  of  free  stream  turbulence  and  skin  friction  was  success¬ 
fully  formulated  by  experiments.  Experimental  results  show  that  Cebeci  &  Bradshaw’s  [59]  mean 
velocity  profile  works  in  the  inner  part.  If  the  free  stream  turbulence  level  is  less  than  12.5%,  rms 
profiles  without  free  stream  turbulence  can  be  applied.  Therefore,  Spalart’s  [52]  inlet  solver  can  be 
safely  utilized  as  long  as  free  stream  turbulence  level  is  lower  than  12.5%. 

But  the  profiles  at  the  outer  part  of  boundary  layer  are  more  complex.  We  take  advantage  of 
Lund  et  al.  ’s  [53]  outer  part  treatment,  but  the  coefficient  is  different, 


awf  =  ^y«c  (HD 

where  S,„/;  is  the  desired  inlet  boundary  layer  thickness,  and  899  is  the  boundary  thickness 
located  at  recycle  plane.  In  [53]  y  was  the  weighting  factor  for  the  outer  part  of  boundary  layer 
associated  with  the  defect  law.  Since  the  defect  law  does  not  work  here,  yi  was  chosen  such  that  the 
inner  part  and  outer  part  was  continuous  in  the  mean  velocity,  i.e. 


j/(y+  =  30) 
Yl  “  Urecy(y+  =  30) 


(112) 


Y2  was  chosen  such  that  rms  velocity  profile  was  continuous  at  the  edge  of  boundary  layer,  i.e. 


TU 

Y2  =  - -T— (113) 

f/"-(y  =  8) 

where  TU  is  the  free  stream  turbulent  level. 

Equation  (110)  represents  our  treatment  of  free  stream  turbulence.  We  need  to  maintain  the  free 
stream  turbulence  level  and  the  turbulent  length  scale.  Since  the  free  stream  does  not  contain  turbu¬ 
lent  production,  small  eddies  decay  faster  than  the  larger  eddies.  Thus,  at  the  different  downstream 
stations,  the  ratio  of  energy  contained  by  small  eddies  and  the  energy  contained  by  large  eddies  are 
different.  When  the  scaled  periodic  boundary  condition  is  utilized,  we  cannot  obtain  a  perfect  inlet 
profile  by  only  recycling  a  single  downstream  plane  because  that  ratio  is  variable.  So  we  recycle 
several  planes,  and  use  those  planes  to  produce  our  inlet  free  stream  profile.  In  Eq.  (110),  weighti 
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Table  3:  Weightings  for  Eq.  (110) 
Station  Location  Weighting 


1  x  =  0.85  1.0/1.75 

2  x  =  1.68  0.5/1.75 

3  x  =  8.05  0.25/1.75 


are  weight  functions  that  monotonically  decrease  in  the  order  of  increasing  downstream  distance, 


^  Weight i  =  1  (114) 

0  <i<m 

and  Y3  is  a  coefficient  such  that  the  TU  is  maintained  in  the  inlet  profile.  This  we  call  a  multi-level 
recycling  method.  The  weighting  function  used  in  this  work  is  given  in  Table  3.  where  3  stations 
were  applied  in  our  calculation.  Such  a  weighting  was  chosen  so  that  the  target  length  scale  would 
be  approached. 


Figure  55:  Comparison  of  mean  velocity  profile  in  a  turbulent  boundary  layer  Re§d  =  2000,  TU  = 
5%  and  Tw  =  Tref,  the  solid  line  is  LES  results,  the  dashed  line  gives  a  DNS  profile  by  Spalart  [52] 
and  the  square  symbols  are  experimental  data  by  DeGraaff  and  Eaton  [60] 

Turbulent  boundary  layers  ranging  from  Red  =  1800  up  to  Red  =  2150  were  calculated  by  a 
dynamic  subgrid  model  proposed  by  Germano  et  al.  [61].  The  simulated  free  stream  turbulent  level 
was  5%.  Figure  55  shows  the  comparison  of  mean  velocity  profile  in  a  turbulent  boundary  layer 
for  Re5d  =  2000  and  TU  =  5%  and  DNS  profile  calculated  by  Spalart  [52]  and  an  experimental 
profile  by  DeGraaff  and  Eaton  [60].  Those  DNS  and  experimental  data  are  all  on  Re§d  =  2000,  but 
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Figure  56:  Comparison  of  streamwise  rms  profiles  in  a  turbulent  boundary  layer  Re§d  =  2000, 
TU  =  5%  and  Tw  =  Tref  The  DNS  results  are  from  Spalart  [52]  and  the  experimental  data  are  from 
DeGraaff  and  Eaton  [60] 

without  free-stream  turbulence.  The  LES  results  in  Fig.  55  exhibit  almost  the  same  main  velocity 
curve  in  the  inner  part  of  boundary  layers  as  the  DNS  and  experimental  data  with  no  free  stream 
turbulence,  but  exhibit  a  larger  difference  in  the  wake  part.  These  phenomena  match  generally 
with  the  experimental  observations  by  Barrett  and  Hollingsworth  [62],  Blair  [55]  and  Hancock  and 
Bradshaw  [56]. 

5.5  Effect  of  Free  Stream  Turbulence  and  Heat  Transfer  on  Turbulent  Boundary 
Layers 

Several  cases  were  calculated  in  order  to  illustrate  the  effect  of  free  stream  turbulence  and  heat 
transfer  on  the  turbulent  boundary  layer.  For  all  cases  the  Mach  number  was  0.06,  and  Re§d  =  2000, 
where  8^  is  the  inflow  displacement  thickness.  The  velocities  in  Fig.  57  were  normalized  by  the  wall 
friction  velocities.  The  figure  shows  the  effect  of  the  heated  wall  on  the  boundary  layer  developing 
with  no  free  stream  turbulence.  With  wall  heating  the  boundary  layer  is  thinner  than  for  the  adiabatic 
case  and  some  deviation  from  the  log  law  is  noted.  The  numerical  results  show  that  heat  transfer 
will  reduce  the  skin  friction,  but  the  free  stream  turbulence  will  increase  it.  Figure  58  illustrates 
the  effect  of  heat  transfer  (without  free  stream  turbulence)  and  free  stream  turbulence  (in  adiabatic 
flow)  to  the  skin  friction  and  Stanton  number.  The  empirical  curve  in  the  Fig.  58  is 
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Figure  57:  Comparison  of  mean  velocity  profiles  in  a  turbulent  boundary  layer  Re =  2000;  the 
solid  and  dotted  lines  are  LES  results,  the  dashed  line  gives  a  DNS  profile  by  Spalart  [52]  and  the 
square  symbols  are  experimental  data  by  DeGraaff  and  Eaton  [60] 

C/  =  0.025/?Cq1/4  (115) 


where  Re%  is  the  Reynolds  number  based  on  momentum  thicknes,  and  Tu  stands  for  turbulence 
level.  Subjected  to  the  free  stream  turbulence  with  5%  turbulence  level,  the  turbulence  length  scale, 
L,  is  495 d  under  the  definition 


(t/2)3/2^'2 

U  dx 


(116) 


where  8</  is  the  inflow  displacement  thickness. 
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Figure  58:  Comparison  of  Cf  and  St  in  a  turbulent  boundary  layer 
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6  Film  Cooling 


This  facet  of  our  work  is  still  in  the  code  development  and  validation  phase.  The  results  to  date 
look  very  promising,  but  additional  simulations  need  to  be  completed  before  final  conclusions  can 
be  drawn.  The  status  of  the  work  is  discussed  below. 

The  effect  of  the  length  of  the  coolant  supply  tube  has  received  considerable  attention  in  the 
literature.  The  tube  L/D  influences  the  tube  discharge  conditions  and  the  resulting  downstream 
structures.  Goldstein  et  al.  [63]  found  an  appreciable  difference  in  effectiveness  between  an  L/D  = 
5.2  and  longer  injection  lengths.  As  the  L/D  increases,  the  discharge  velocity  profiles  become  more 
uniform.  For  a  long  tube,  say  L/D  >  6,  the  discharge  condition  becomes  nearly  independent  of  tube 
entrance  conditions. 

Burd  et  al.  [64]  compared  the  film  cooling  effectiveness  of  an  L/D  =  7.0  configuration  with  one 
having  L/D  —  2.3.  Their  experiments  showed  that  the  short-hole  injection  flow  penetrates  farther 
from  the  wall  and  influences  a  greater  extent  of  the  region  downstream  from  the  hole  under  low  free 
stream  turbulence,  but  no  significant  differences  in  the  normalized  mean  velocities  were  noted  with 
high  free  stream  turbulence. 

Ligrani  et  al.  [65]  examined  the  film  cooling  effectiveness  of  a  single  row  of  film-cooling  holes 
with  and  without  compound  angle  orientations.  They  evaluated  blowing  ratios  of  0.5,  1.0  and  1.5. 

A  long-tube  configuration  with  L/D  —  8  was  selected  for  the  first  LES  film  cooling  simulation. 
The  blowing  ratio  was  0.5  and  the  density  ratio,  2.0.  The  wall  was  adiabatic.  The  arrangement  of 
the  computational  domain  is  shown  in  Fig.  59.  The  ratio  of  the  diameter  of  pipe  to  the  width  of 


Figure  59:  The  numerical  configuration  for  film  cooling 


boundary  layer  namely  D/Z  was  1  /5.  The  ratio  of  the  diameter  of  pipe  to  the  boundary  thickness, 
namely  D/8,  was  0.44.  The  Reynolds  number  of  the  boundary  layer  ranged  from  1800  up  to  2700 
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based  on  the  displacement  thickness.  The  distance  from  inlet  to  the  center  of  the  coolant  discharge 
hole  was  56.58^,  where  S<*  is  displacement  thickness,  and  the  distance  from  the  hole  center  to  the 
outflow  boundary  was  72. 58,/,  which  which  corresponds  to  22.656 D.  The  Reynolds  number  of  the 
boundary  layer  at  the  discharge  hold  center  was  2000.  The  Reynolds  number  for  the  pipe  flow  was 
9357.1.  The  injection  angle  was  35°.  The  film  cooling  effectiveness  was  calculated  according  to 


ri  = 


(117) 


The  grid  was  non-orthogonal  in  the  tube  and  near  the  injection  hole.  Conservation  was  main¬ 
tained  by  employing  control  volume  concepts  similar  those  employed  with  unstructured  meshes. 
Figure  60  provides  a  downward  view  of  the  numerical  mesh  employed  for  the  external  boundary 
layer  near  the  injection  hole.  Figure  61  shows  the  mesh  in  more  detail. 


Figure  60:  The  computational  mesh  near  the  injection  hole 

Sixteen  processors  were  used  to  calculate  this  case.  The  LU-SGS  scheme  and  a  dynamic  SGS 
model  were  utilized.  The  numerical  grid  for  the  external  turbulent  boundary  layer  was  591  x 
85  x  101  corresponding  to  the  streamwise,  normal  and  spanwise  directions,  respectively.  The  nu¬ 
merical  domain  was  126  x  22  x  16  in  units  of  displacement  thickness.  The  numerical  grid  was 
241  x  41  x  101  corresponding  to  streamwise,  normal,  and  spanwise  directions  for  the  coolant  supply 
pipe, respectively.  The  total  number  of  control  volumes  used  in  the  simulation  exceeded  6  million. 
In  units  of  the  external  boundary  layer  initial  displacement  thickness,  the  pipe  radius  was  1.6  and 
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Figure  6 1 :  Expanded  view  of  computational  mesh  near  the  leading  edge  of  the  injection  hole 


the  length,  16. 
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Figure  62:  The  downward  view  of  instantaneous  temperature  contours  at  =  14;  the  units  of  the 
axes  are  displacement  thickness,  and  the  center  of  the  hole  is  located  at  x  =  56.5;  Y  represents  the 
spanwise  direction 


Figure  62  shows  the  evolution  of  the  coolant  in  the  buffer  zone  of  the  turbulent  boundary  layer. 
Due  to  the  strong  turbulence  in  this  region,  the  coolant  mixes  with  the  external  flow  very  rapidly. 
Such  mixing  accelerates  the  energy  transfer  between  the  coolant  and  the  surrounding  flow. 

Figure  63  shows  the  side  view  of  the  instantaneous  temperature  contours.  This  figure  clearly 
shows  the  evolution  of  the  coolant  downstream.  Although  not  discemable  in  the  temperature  field, 
the  flow  reverses  immediately  downstream  of  the  hole.  Such  a  reversed  flow  can  be  visualized  as  a 
vortex.  The  coolant  reattaches  to  the  wall  just  behind  this  vortex.  The  reattachment  accelerated  the 
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Figure  63:  The  side  view  of  instantaneous  temperature  contours  at  the  centerline  of  the  numerical 
domain,  the  units  of  the  axes  are  displacement  thickness,  and  the  center  of  the  hole  is  located  at 


x  =  56.5 


mixing. 

Figure  64  shows  the  distribution  of  adiabatic  effectiveness  on  the  flat  plate.  In  this  case,  the 
Reynolds  number  of  the  external  turbulent  boundary  layer  ranges  from  1800  to  2300  based  on  dis¬ 
placement  thickness,  the  density  ratio  is  2.0,  and  the  blowing  ratio  is  0.5.  The  free  stream  turbulence 
level  is  zero  in  this  case.  Figure  64  shows  that  the  distribution  of  adiabatic  effectiveness  on  the  wall 
is  nearly  symmetric  about  the  center  line  of  the  hole.  The  magnitude  decays  in  the  streamwise 
direction  as  the  cooled  area  increases  in  width. 

Figure  65  shows  the  distribution  of  adiabatic  effectiveness  on  the  wall  at  the  station  X/D  — 
8.75  downstream  of  the  hole.  This  plot  shows  a  symmetric  distribution  which  agrees  reasonably 
with  experimental  results  [66],  [65]  and  [64].  Figure  66  shows  the  corresponding  contours  of  the 
instantaneous  temperature  field  at  this  station. 

Figure  67  is  the  plot  of  adiabatic  effectiveness  along  the  center  line  downstream  of  the  hole. 
Very  near  the  hole,  the  adiabatic  effectiveness  decreases  due  to  a  region  of  reversed  flow.  The 
local  minimum  was  located  at  X/D  =  0.9  from  the  center  of  the  hole.  The  peak  value  occurred 
at  X/D  =  3.2  from  the  center  of  the  hole.  The  increase  in  adiabatic  effectiveness  results  from 
the  reattachment  of  the  flow,  which  increases  the  heat  transfer  at  the  wall.  The  peak  value  of  the 
adiabatic  effectiveness  for  this  case  with  a  density  ratio  of  2.0  and  a  velocity  ratio  of  0.5  was  0.78. 
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Figure  64:  The  distribution  of  adiabatic  effectiveness  on  the  flat  plate  for  the  case  of  density  ratio 
DR  =  2.0,  blowing  ratio  VR  =  0.5,  free  stream  turbulence  level  Tu  =  0.0;  X  is  the  streamwise 
direction  and  Y  is  the  spanwise  direction 


The  LES  results  for  film  cooling  look  very  promising.  However,  much  more  remains  to  be 
done  in  order  to  validate  the  procedure  and  apply  it  to  study  important  issues  related  to  physics  and 
design.  The  remaining  task  of  highest  priority  is  to  make  comparisons  with  experimental  results. 
A  case  with  a  density  ratio  of  two  was  chosen  in  order  to  emphasize  the  power  of  the  compressible 
LES  procedure.  Unfortunately,  we  know  of  no  experimental  results  for  such  a  large  density  ratio. 

As  funding  permits,  film  cooling  work  that  will  include  the  use  of  shorter  supply  tubes  and  a 
supply  plenum  in  the  computational  domain  will  be  carried  out.  Also  of  considerable  interest  is  the 
inclusion  of  significant  free  stream  turbulence.  The  use  of  pulsed  (unsteady)  coolant  supply  is  also 
of  interest. 
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Figure  65:  Plot  of  adiabatic  effectiveness  at  X/D  =  8.75  downstream  of  the  hole  for  density  ratio 
DR  =  2.0,  blowing  ratio  VR  =  0.5  and  free  stream  turbulence  level  Tu  =  0.0 
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Figure  66:  Upstream  view  of  instantaneous  temperature  contours  at  X=  84  (X/D  =  8.75);  the  units 
of  the  axes  are  displacement  thickness 
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Figure  67:  Adiabatic  effectiveness  along  the  center  line  of  the  flat  plate  for  density  ratio  DR  =  2.0, 
blowing  ratio  VR  —  0.5  and  free  stream  turbulence  level  Tu  =  0.0 
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7  Conclusions 


Large  eddy  simulation  has  been  shown  to  be  a  very  useful  tool  to  aid  understanding  of  the 
complex  flows  characteristic  of  turbine  blade  cooling  systems.  At  this  point  the  LES  technology 
is  limited  to  relatively  low  Reynolds  numbers  and  simple  geometry,  but  effects  due  to  rotation  and 
heat  transfer  can  be  studied. 

The  major  conclusions  from  this  study  are  given  below. 

•  Although  still  somewhat  limited  to  low  Reynolds  numbers  and  relatively  simple  geometries, 
large  eddy  simulation  was  observed  to  be  a  valuable  tool  for  examining  the  behavior  of  tur¬ 
bulent  flows  with  rotation  or  rib-roughness  that  are  difficult  to  compute  accurately  by  other 
approaches. 

•  In  general,  system  rotation  was  found  to  suppress  turbulent  velocity  fluctuations  and  shear 
stresses  near  the  leading  side  of  a  rotating  channel  and  enhance  them  near  the  trailing  side. 
The  ratio  of  the  skin-friction  coefficient  and  Nusselt  numbers  on  the  two  sides  ranged  between 
a  factor  of  2  and  3.  Heat  transfer  appeared  to  be  influenced  somewhat  more  strongly  than  skin 
friction.  The  maximum  mean  streamwise  velocity  was  skewed  toward  the  leading  side  of  the 
rotating  channel. 

•  In  the  rotating  ribbed  channel,  rotation  caused  the  length  of  the  separated  flow  region  behind 
the  ribs  to  lengthen  on  the  leading  side  and  shrink  on  the  trailing  side.  This  appears  consis¬ 
tent  with  the  observed  trend  in  shear  stresses;  i.e.,  increased  stresses  resulted  in  more  rapid 
reattachment. 

•  For  a  rotating  square  duct,  a  distinct  secondary  flow  pattern  was  observed  in  the  simulations 
in  general  agreement  with  experimental  observations.  Flow  near  the  sidewalls  was  driven  by 
the  pressure  gradient  from  the  trailing  side  to  the  leading  side  and  a  less  intense  flow  was 
observed  in  the  central  region  from  the  leading  side  to  the  trailing  side.  The  latter  migration 
resulted  in  the  maximum  in  the  mean  streamwise  velocity  distribution  being  shifted  toward 
the  trailing  side,  opposite  from  the  skewing  observed  in  channel  flow. 

•  Centrifugal  buoyancy  terms  were  included  in  the  governing  equations  for  the  rotating  square 
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duct  flows  with  heat  transfer.  The  centrifugal  effect  was  observed  to  significantly  alter  the 
flow  pattern.  When  Gr  <  0,  the  buoyancy  force  was  observed  to  cause  flow  separation  on  the 
leading  wall. 

•  Work  was  completed  toward  an  efficient  fractional  step  method  for  solving  the  compressible 
Navier-Stokes  equations  that  has  some  novel  features.  The  scheme  was  used  in  the  work  on 
external  turbulent  boundary  layers  that  was  necessary  for  the  effort  on  film  cooling. 

•  A  procedure  was  developed  to  establish  inflow  conditions  for  LES  of  spatially  developing 
turbulent  boundary  layers  that  included  the  influence  of  freestream  turbulence.  This  also  was 
part  of  the  film  cooling  effort. 

•  A  demonstration  simulation  for  an  idealized  film  cooling  configuration  was  completed.  The 
essential  components  for  film  cooling  simulations  are  in  place  and  the  preliminary  results 
appear  reasonable  and  promising.  Work  is  underway  to  make  quantitative  comparisons  with 
experimental  data  to  validate  the  approach  and  confirm  the  level  of  accuracy. 
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